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Abstract 

Convection-diffusion equations provide the basis for describing heat and mass 
transfer phenomena as well as processes of continuum mechanics. To handle 
flows in porous media, the fundamental issue is to model correctly the convective 
transport of individual phases. Moreover, for compressible media, the pressure 
equation itself is just a time-dependent convection-diffusion equation. 

For different problems, a convection-diffusion equation may be be written in 
various forms. The most popular formulation of convective transport employs the 
divergent (conservative) form. In some cases, the nondivergent (characteristic) 
form seems to be preferable. The so-called skew-symmetric form of convective 
transport operators that is the half-sum of the operators in the divergent and non- 
divergent forms is of great interest in some applications. 

Here we discuss the basic classes of discretization in space: finite difference 
schemes on rectangular grids, approximations on general polyhedra (the finite 
volume method), and finite element procedures. The key properties of discrete 
operators are studied for convective and diffusive transport. We emphasize the 
problems of constructing approximations for convection and diffusion operators 
that satisfy the maximum principle at the discrete level — they are called mono- 
tone approximations. 

Two- and three-level schemes are investigated for transient problems. Uncon- 
ditionally stable explicit-implicit schemes are developed for convection-diffusion 
problems. Stability conditions are obtained both in finite-dimensional Hilbert 
spaces and in Banach spaces depending on the form in which the convection- 
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diffusion equation is written. 
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1. Convection-diffusion problems 

Convection-diffusion problems are governed by typical mathematical models, 
which are common in fluid and gas dynamics. Heat and mass transfer is con- 
ducted not only via diffusion, but appears due to motion of a medium, too. Here 
we present typical examples of model convection-diffusion problems, which use 
various forms for the terms describing convective transport. 

1.1. Basic problems of continuum mechanics 

Principal features of physical and chemical processes in fluid mechanics ifTil 
3] result from motion of a medium due to various forces. Heat and mass transfer 
phenomena in a moving medium may be treated as the simplest examples of these 
peculiarities. 

Let v(x,t) be the velocity of a liquid at a point x and at a time moment t, 
whereas p is its density. The thermal state of the liquid is governed by the equation 
of heat conduction 



where T stands for the temperature, Cp is the specific heat capacity at a constant 
pressure, k denotes the thermal conductivity of the liquid and q describes the 
intensity of volumetric heat sources. 

The temperature at a given spatial point is governed not only by conduction 
(diffusion) of heat, but also by motion (convection) of fluid volumes. 

The second typical example is the diffusion equation for a multicomponent 
mixture [jsi @]- We assume that a liquid is heterogeneous, more exactly, it is a 
mixture of two components. In this case, the mixture composition may be de- 
scribed by the concentration c associated with only one component. The corre- 
sponding equation for the concentration (neglecting the diffusion flux caused by 
the temperature gradient) has the form 




(1) 



d{pc) 



+ grad(vpc) = div(pA; grade), 



(2) 



dt 
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where k denotes the diffusivity and p is treated as the total density of the liquid. 
The equation (O may be rewritten as 

— — h divivm) + div I m— gradp ) = div(A;gradm), (3) 
dt \ P J 

where m = pc is the mass of one of the components in a volume unit. The 
equation ^ may be reduced to the form 

— — + div(Dm) = div(A;gradm), (4) 

where the expression 

V = V -\ — grad p 

P 

describes the effective convective transport. 
Using the continuity equation 

|^ + divM = 0, (5) 

we obtain 

+ div(vpc) = p[— + {v grad)c 

Therefore, under the natural assumptions on the positiveness of p, from (O, we 
arrive at the equation 

dc k 

— + {v ■ grad)c grad p ■ grad c = div(fc grad c). (6) 

at p 

Similarly to dU, rewrite ^ as 

dc 

— + (?;■ grad)c = div( A; grade), (7) 
ot 

where now 

V = V grad p. 

P 

Thus, we come to the equation for the concentration, where convective transport 
has the nondivergent form, as it takes place in the heat equation ([T]). In equation 
dH) as well as in the continuity equation dS]), convective transport is written in the 
divergent form. 
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More complete models of heat and mass transfer include also an equation that 
describes the motion of the medium itself and determines, in particular, the veloc- 
ity V. For simplicity, we restrict ourselves to the Navier-Stokes equation for an 
incompressible (p = const) homogeneous medium. In this case, the momentum 
equation seems like this: 

fdv \ 
p I — + {v ■ gTa.d)v = — gradp + r] div grad v, (8) 



\dt 

whereas the continuity equation ^ is reduced to 

div v = 0. (9) 

Here p denotes the pressure and t] = const stands for the dynamic viscosity of the 
fluid. 

The equations ([8]) may be treated as the equations of convective and diffusive 
(due to the viscosity) transport of each individual component of the velocity v. In 
this situation, in order to evaluate the pressure p, it is necessary to involve equation 

m- 

If we eliminate the pressure from equation dS]), then, for the vorticity w = 
rot V, we obtain the equation 

(dw \ 
+ {v ■ grad)it> — {w ■ grad)f j = rj div grad w. 

It is easy to see that the dynamics of the vorticity for an incompressible fluid is 
determined by a specific convective and diffusive transport. 

More sophisticated models that include convective and diffusive transport as 
the most important element are used in modeling compressible flows. It should be 
noted that convective-diffusive transport is essential for predictions of various gas 
and fluid flows. In particular, environmental problems are of great importance: 
pollutants spreading in the atmosphere and water basins, contaminants transport 
in groundwaters and so on. 

1.2. Various forms of the hydrodynamics equations 

In theoretical studying applied problems, the conservative form of the hydro- 
dynamics equations is in common use. In this case, the equations have the di- 
vergent form and express directly the corresponding laws of conservation (for the 
mass, momentum and energy). On the other hand, we should pay attention to 
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the nondivergent (characteristic) form of the hydrodynamics equations, which is 
connected with the representation that is derived via differentiating the convective 
transport terms. The paper [i38i1 presents a new form of the hydrodynamics equa- 
tions that is characterized by writing the convective terms in the skew- symmetric 
form. New quantities — the so-called SD-variables (Square root from Density) 
that are based on using not the density itself but the square root from the density 
— are used as unknown variables. Physical and mathematical arguments in favor 
of introducing this form of the hydrodynamics equations are discussed below. 

The system of hydrodynamics equations includes, first of all, the scalar equa- 
tion of continuity and the vector equation of momentum. In more common cases, 
there may be several motion equations as well as continuity equations — we speak 
of models for multicomponent media. Furthermore, the system of equations may 
be supplemented with an energy equation. Usually, the following scalar equation 
of convection-diffusion serves as the basic equation in continuum mechanics and 
heat and mass transfer (see, e.g., [20,|43|]), i.e., 

+ diY(gvip) = div(Dgrad(^), (10) 

at 

where is a desired scalar function and D denotes the diffusivity. This equa- 
tion is written in the conservative (divergent) form. Concerning to equation (flOl) . 
a number of problems are discussed in the literature, such as the construction 
of discretization in space and in time, the investigation of convergence of the ap- 
proximate solution to the exact one for the corresponding boundary value problem 

HQ. 

The main peculiarities of the system of fluid dynamic equations become evi- 
dent in studying the system of two scalar equations that includes not only equation(fTOl). 
but also the continuity equation 

|^ + div(H = o. (11) 

Just this system of equations (fTOl) . ([TT]) is said to be the basic system of scalar 
hydrodynamic equations. 

In investigation of transport phenomena in continuum mechanics, the primary 
features of transport of scalar variables are represented in equation (flOl) . Concern- 
ing vector fields, the coordinatewise representation may be unsuitable. Thus, it 
seems reasonable to supplement the system of equations (fTOl ), (fTTI) with the vector 
convection-diffusion equation 

d{Qu) 



dt 



+ div(f)i> ® u) = div(D grad w), (12) 
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where u is the desired vector function. This system of equations (fT0l)-(fT2l) is 
called the basic system of hydrodynamics equations. 
Taking into account 

dw{gv(f) = (pdw{gv) + gv ■ grady? 

and the continuity equation (fTTI) . we get 

dip 

g— — \- gv ■ gradip = div{D gra.d{p). (13) 



Similarly, we can rewrite equation (|T2l) as 



g— — h ^''y • grad n = div(D grad n). (14) 



The equations (1131) . (1141) are written in the nonconservative (nondivergent) form. 
It should be noted that the continuity equation (fTTI) cannot be written in the non- 
conservative form. Therefore, the basic system of hydrodynamics equations may 
be written in the conservative (fT0l) - (fT2)) or in the partially nonconservative form 
(fTTI) . (fT3l) . (fT4l) . Only for an incompressible medium, where equation (fTTI) takes 
the form 

div V = 0, 

it is possible to speak about the nonconservative form of the equations. 

Let us write the operator of convective transport in the skew- symmetric form 
HH as 

Ce = ^ diY{v9) + ■ grad 9, (15) 

i.e., as the half-sum of the operators of convective transport in divergent (conser- 
vative) and nondivergent (nonconservative) forms. 

In the basic system of fluid dynamics equations (fT0l)-(fT2l). instead of g, ip, u, 
we introduce new unknown variables: 

s = {gY/', (={qY^'v, w = {gy^'u. (16) 

The main peculiarity of these unknowns consists in using the square root from 
the density s = (gY^"^ instead the density g itself. That is why we speak of SD- 
variables (Square root from Density). 

For the new unknowns, the system of equations (fT0|)- (fT2l) may be rewritten in 
the following way: 

ds 1 1 

— + - div(vs) + -V ■ grad s = 0, (17) 
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^ + 1 div(i;C) + \v ■ grade = \ div f/^grad (^]], (18) 



+ - div(v ® w) -\- -17 ■ grad t/; = - div (^D grad ^ — j j . (19) 



In this case, all three equations involve the convective terms that are written in the 
skew- symmetric form. 

As a typical example of using the new variables, we study the Navier-Stokes 
equations for a viscous compressible medium, which express the conservations 
laws for the mass, momentum, and energy. The continuity equation has the form 
(fTTI) . Usually, the momentum equation is written in the conservative form 

— — h (iiv{Qv ® v) = div N — gradp. (20) 

at 

Here 

2 

N = — /idiv V E + 2/iS, 

3 

As for S, the coordinatewise representation seems like this: 

S.. = i + ^ 

2 \dxj dxi 

Now introduce the energy equation 

d{Qe) 



dt 



+ div(f)'ue) = div(fcgradT) + pdhrv + N : gradt;. (21) 



The term N : grad v describes the heat dissipation due to the fluid viscosity and 
N : grad v is the scalar product of tensors: 

N : grad v= N^x^ + Nxy^ + Nx,^ 
ox Oy Oz 

+N ^ + iV ^ + iV ^ 

'■'-^yx o ' -'■"yy r% ~ -'■"yz ^ 

OX oy oz 

dvz dvy ,^ 

ax ay az 

Let us introduce the following new unknown variables: 

s = {ef'\ w={gY/'v, C = {Qf"e. ill) 
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For the variables (l22l) . the system of the Navier-Stokes equations (fTTl) . (|20l ). (|2T]) 
has the following form: 

dw 1/ /W \ w \ l,., 1 

+ - ( div ( — (8) w ) H ■ grad w j = - div N grad p, (24) 



dt 2 V V s / s 
d( 1 / \ w 



dt 



1 p /w\ 1 /w 

- div(A;gradT) H — div — ) H — N : grad — 
s s \ s / s \ s 



(25) 



The system of equations (|23l) - (|25l) needs to be supplemented with some equation 
of state. It should be highlighted that using the variables ((22)) . the convective terms 
are written in the skew- symmetric form. 

1.3. The pressure problem for multiphase flows in porous media 

The system of governing equations for multiphase flows includes [21, 1] the 



continuity equation for each phase, where a = 1, 2, . . . , m is the phase index. The 
mass conservation law for each individual phase is expressed by the following 
equation: 

— h dw{baUa) = -baQa, a = 1, 2, . . . , m. (26) 

Here stands for the porosity, ha is the phase density, Sa denotes the phase satu- 
ration, Ua is the velocity, and Qa describe the volumetric mass sources. 

For simplicity, we neglect the capillary and gravitational forces. In this sim- 
plest case, the equation of fluid motion in porous media has the form of Darcy's 
law, where the velocity is directly determined by the common pressure: 

k 

Ua = — - k ■ gradp, a = 1,2, . . . ,m. (27) 

In (ITT] ), k is the absolute permeability (in general, a symmetric second-rank ten- 
sor), ka denote the relative permeabilities, /Xq, stands for the phase viscosity, and 
p is the pressure. 

The unknown variables in the system of equations (|26l) . (1271) are the phase 
saturations Sa, a = 1,2, ... ,Tn and the pressure (m + 1 unknowns in all). In the 
simplest case, the coefficients in equations (|26l ). (ITTI) are defined as some relations 

(j) = (f){p), ba = baip), Qa = QaiSa), k^ = kaiS^), fJ^a = COHSt . 
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The summation of saturations over all phases yields 

m 

Y,Sa = l. (28) 

a=l 

Substituting (|271) in (|26l) and taking into account (|28] ). we obtain a system of m + 1 
equations for m + 1 unknowns. 

The system of equations (I26l)-(l28l) provides the basis for the description of 
multiphase flows in porous media. We have no separate equation for the pressure 
in this system. The equations (|26l) may be treated as the transport equation for 
each individual phase, whereas the algebraic relation (|28l) may be considered as 
the equation for the pressure. 

Let us consider more convenient forms for the system (I26l)-(l28l). which lead 
to the typical problems of mathematical physics for the pressure. It should be 
noted that such equivalent formulations do exist only at the differential level. At 
the discrete level, such equivalence of formulations is not valid even for linear 
problems. Thus, a proper choice of the original form of the governing equations 
is essential for calculations. 

The most natural way to derive the equation for the pressure is the following. 
Dividing each equation (l26l) by 6q, > and adding them together, we get 

Sa d{(j)ba)\ dp 1 fbaka. .\ 1 .^Q. 

(^E ^— j = E ^ d,v k . gxadpj - - g (29) 

Under the natural assumption for compressible fluids that 

; > 0, a = l,2,...,m, 

dp 

equation (|29l ) for the pressure is the standard parabolic equation of second order. 
In particular, the maximum principle holds for its solutions [10 ]. 

In accordance with (|29l) . we solve in the boundary value problem for the 
equation 

du 

— + '^aa{x)CaU = f{x,t), (30) 

where aa{x) > Qa, Qa > 0, a = 1,2, ...,m, and the elliptic operators are 
defined by 
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under the standard assumptions < < fca < Kq,. 

In some cases (incompressible media), it is reasonable to consider the steady- 
state problem. The boundary value problem is formulated for the equation 

m 

Cu = f{x), C = ^aaix)Ca, (32) 

a=l 

which is supplemented by the boundary conditions. 
From (|3T|) and (|321 ). we have the representation 

m 

C = ^Ca, Ca=Va+Ca, a = 1,2, ...,171, (33) 

a=l 

where 

VaU = — div{da{x) grad-u), (34) 
CoM = Wagradu. (35) 
The effective diffusivity and convection velocity for the individual phase a are 

Then the pressure operator takes the form of the convection-diffusion operator 
with convective term in the nondivergent form. 

2. Time-dependent problems of convection-diffusion 

Convection-diffusion equations provide important examples of second-order 
parabolic equations. In particular, they are considered as the basic equations for 
modeling continuum mechanics phenomena. Some aspects of numerical solving 
time-dependent problems of convection-diffusion are discussed here. In these 
equations, convective terms are formulated in the divergent, nondivergent, and 
skew- symmetric forms. Some essential results are presented for a model initial- 
boundary value problem with Dirichlet boundary conditions for the differential 
equation of convection-diffusion. These results will serve us as a check point in 
developing difference schemes. Discrete operators of diffusion and convection are 
constructed and analyzed with respect to their primary properties. 
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2.1. Differential problems 

Time-dependent problems of convection-diffusion are treated as evolutionary 
operator equations in the corresponding spaces. To investigate them, we start 
with a study on properties of differential operators describing convective and dif- 
fusive transport. As the basic problem, we consider a time-dependent problem 
of convection-diffusion with Dirichlet boundary conditions in a rectangle. The 
convective terms are written in various forms. We distinguish a class of model 
time-dependent problems of convection-diffusion with a constant coefficient of 
diffusive transport (it is independent of time but depends on spatial coordinates). 
As for coefficients of convective transport, in applied problems, they are variable 
both in space and in time. 

In a rectangle 

n = {x\x = (Xi, X2) , < Xa < la, a — 1, 2}, 

we study the time-dependent convection-diffusion equation with the convective 
terms written in the nondivergent form: 

du v-^ . V du 



dt dxa 

(36) 

considered under the standard assumptions ki < k {x) < k^, ki > 0, T > 0. 
This equation is supplemented with homogeneous Dirichlet boundary conditions 

u{x, t) = 0, xedVt, <t<T. (37) 

For the unique solvability of the unsteady problem, it is supplemented with the 
initial condition 

u{x,0) = u^{x), xen. (38) 

The second example is the time-dependent equation of convection-diffusion 
with the convective transport written in the divergence form: 

du 9 , , , , 

-1 (39) 
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And finally, the primary object of our investigation is the convection-diffusion 
equation with the convective terms written in the skew-symmetric form: 



(40) 

a=l V / 

We consider a set of functions u{x) that satisfy the boundary condition (l37l) . 
The transient convection-diffusion problem may be formulated as the operator- 
differential equation 

fjii 

— + Au = f{t), A = C + V. (41) 
Here V is the diffusive transport operator that is defined by the expression 



According to (|36l) . (1391) . (1401) . the convective transport operator is written in 
distinct forms. For the convective transport operator in the nondivergent form, 
from (l36l) . we set C = Ci, where 

^ du 



C^u = ^v^{x,t)-^. (43) 

Similarly, from (|39| ). we have C = C2, where now 

2 



C2M = ^^K(a;,t)M). (44) 

Taking into account (|40l ), the convective-transport operator in the skew-symmetric 
form is 

C = Co = -(C1 + C2), 

and 

Now we highlight the basic properties of the above-mentioned operators of diffu- 
sive and convective transport. 
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2.2. Properties of the operators of diffusive and convective transport 

The solution of a discrete problem should inherit the basic properties of the 
corresponding differential problem. This can be achieved, in particular, if the grid 
operators have the same primary properties as the differential ones. 

As usually, let = iv2(f2) be a Hilbert space for arbitrary functions u(x) 
and w{x) equal zero on dil. The diffusive transport operator defined by (|42|) is 
self-adjoint in T-L on the set of functions satisfying the boundary conditions (|37l) : 

V = V*. (46) 

Note also that the diffusive transport operator under consideration at the above- 
mentioned restrictions is positive definite, i.e., the estimate 

V > KiXoS, (47) 

is valid, where S denotes the identity operator and Aq > is the minimal eigen- 
value of the Laplace operator with the Dirichlet boundary conditions. For the 
rectangle il, we have 

The estimate (|47|) follows from 



(Vu,u) > Ki(Vu,'Vu) > KiXo(u,u). 

We now consider the convective transport operator in various formulations 
(see (|43l) . (|44l) . and (l45l)). Taking into account the homogeneous boundary condi- 
tions (l37l) . we have 



(CiU,w)=y / VaT^ Wdx = —y / — (VaW)udx = —(U,C2W). 



Thus, we see that the convective transport operators in the divergent and nondi- 
vergent forms are the adjoints of each other (with a precision of the sign): 

C* = -C2. (48) 

In view of (l48l) . the convective transport operator in the skew- symmetric form 
(|45l) is skew-symmetric ((CqM, u) = 0): 

Co = -C*. (49) 
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Under the condition of incompressibility 

2 o 

div v = y = 0, xen, (50) 

a=l 

the convective transport operator in the nondivergent fl3l) and divergent (|44|) forms 
are also skew- symmetric. In constructing discrete approximations of the convec- 
tive transport operators, the principal moment is that the skew- symmetric property 
of the operator Cq is valid for any Va{x,t), a = 1,2 including the compressible 
case. 

It seems useful to give the upper bound for the convective transport operator. 
For (03]), (gH), we have 

X ^ f T f 

{Ciu, u) = —{C2U, u) = - / Va-^ — dx = — - / div v dx 

and therefore 

|(C„m,m)| < ^11 div'y||c(f7) IImT, a = 1,2, (51) 

where 

||u;||c(f7) = max \w{x)\. 

Thus, for the convective transport operators defined in accordance with (|43] ), (|44l) 
(C = Cq, a = 1, 2), we obtain 

\{Cu,u)\<Mi\\u\\\ (52) 

where a constant Mi depends only on div v and, in accordance with (|5T|) , it fol- 
lows that ^ 

Mi = -\\diYv\\c{n)- (53) 

In addition, we present the estimates of subordination of the convective trans- 
port operator to the diffusive transport operator: 

WCuW'^ <M2{Vu,u), (54) 

with a constant M.2 depending on the velocity. 
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For the nondivergent operator of convection (|43l) . we have 

\\cM?= I l> dx<2y]ivi[i^] dx 



du ^ ^ 



<2max{K||c(.)}-$:yA:^^j dx 

(X — 1 (-^ 



2 

< — max{||t;^||c(Q)} {Vu,u) 

/tl a 



i.e., in the inequality (1541) . for C = Ci, we can assume 

2 

A^2 = — max {II 17^11 . 
Similarly, for C = C2 (see dS])), we obtain 

||C2ii|P = J i + div vu j da; 



f7 ^"=1 

2 

\2„,2, 



dx + 2 J (divv) u dx. 
n ' o 

Taking into account the Friedrichs inequality 

2 ^ / o \ 2 



where the constant A^o = ^, we obtain at C = C2 the estimatr (|54l) with 



^2 = ^ (^2max{||i;^||c(n)} +A^o|| div 'i;||^(f2)) . 



2 

In a similar way, for the case C = Co, we have 



I.e., 



M2 = - f3max{||t;^||c(m} + A^oll divw| 



K 

15 
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The above estimates (|52l) . (|54l) serve as a reference point in studying discrete 
analogs of the convective transport operator. 

Summarizing the above properties of the convective transport operator, we 
obtain the following statement. 

Theorem 1. The convective transport operators have the following properties: 

• the convective transport operators in the divergent and nondivergent forms 
are adjoint to each other up to the sign — the equality 

• the convective transport operator in the skew-symmetric form is skew-symmetric 
— the equalitY^49\); 

• the convective transport operators in the divergent and nondivergent forms 
are bounded — the a priori estimates ( |52]) and ([53\) ; 

• the convective transport operators are subordinated to the diffusion oper- 
ator - the estimate 041 ) with the constant M.2, defined according to 051) . 



It seems reasonable to construct difference operator of convective and diffu- 
sive transport in such a way that they do have the same properties. 

2.3. A priori estimates 

We restrict ourselves to elementary a priori estimates for the time-dependent 
equation (14T1) supplemented with the initial condition 



They are based on the above-established properties (see Theorem [T]) of the opera- 
tors of diffusive and convective transport. 

Theorem 2. For the solution of the problem f l471) . 091) . under the conditions f l47l) . 
021 ). 0?1) . the following a priori estimate holds: 



(112), (HI. 



m(0) = 



(59) 



u{t)\\^ < exp(2Mii)||Mo 




(60) 
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\u{t)\\ < exp{^M2t)\\uo\\ 
t 

+ j eM\M2{t-e))\\fmde, 



|Vn(t)f <^exp(>l2t)||V«oir 

t 

+ - [ exp{M2{t-emf{e)fde, 



K 







where ki < k{x) < ^2, x E QU dQ. 
Proof. Multiplying equation (|4T]) scalarly by u{t), we get 
1 d 



Taking into account (1521) and the inequality 

1 



from (l63l). it follows that 



(61) 



(62) 



+ iVu,u) = -{Cu,u) + if,u). (63) 



{f,u)<iVu,u) + ^\\f\\U, 



j^\\ur<2M4ur + l\\f\\l... 
Using Gronwall's lemma, we obtain from this inequality the required estimate 
From (|54l) . we have 

I - {Cu,u)\ < < (Vu,u) + ^M2\\uf. 



This allows to obtain from (1631) the inequality 

Ju\\ < -^M2\\u\ 
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which immediately implies the estimate (|6T1) . 

It remains to derive the estimate (|62l) . To do this, multiply equation PT]) 
scalarly by du / dt and obtain 



du 



dt 



du 



For the right-hand side, we have 

f ^ du\ ( , du\ 

In view of (l54l) . we get the inequality 



du 



dt 



, du 



Incur +1 



By 



< M2iVu,u) + 



kiIIVmIP < (T>u,u) < K2\\Vu\ 



(64) 



from (1641) . we obtain the desired estimate (|62l) . □ 
The resulting estimates (l60l)- (l62l) provide the continuity of the solution of 
((4T)) . (l59l) with respect to the initial data and the right-hand side. In these esti- 
mates, the essential issue is that for the solution norm of the problem with the 
homogeneous right-hand side, it is allowed an exponential growth with a growth 
increment that depends on the constants Aii, M.2- It is necessary to allow such 
a behavior for the solution at the discrete level. Thus, we need to introduce the 
concept of f?- stability for the corresponding difference schemes. 

2.4. The maximum principle and a priori estimates 

Considering boundary value problems both for parabolic equations of the sec- 
ond order in space and for second-order elliptic equations, special attention is 
paid to the maximum principle uM . The corresponding statement is formulated 
as follows. 

Theorem 3. Assume that in the Cauchy problem d47]) . d59l) the right-hand side 
f{x,t) > {f{x,t) < 0) and the initial data uo{x) > {uq{x) < 0), then 
u{x, 0) > {u{x,0) < 0) for all x e and t > 0. 
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Note that it is possible to use the maximum principle in a stronger form that 
employs weak inequalities for the right-hand side, i.e., the non-negativity of the 
solution takes place under the condition of the non-negativity of the right-hand 
side and the initial data. 

Here are some a priori estimates for the convection-diffusion problem (|4T]) . 
(|59l ), which are derived via the maximum principle. In the above-considered time- 
dependent problems with Dirichlet boundary conditions, we can easily construct 
a majorant function — a wide range of the appropriate estimates is given, e.g., 
in the book [13]. We also give an estimate for the convection-diffusion equation 
with convective terms in the divergence form — the estimate in Ci^l). 

Theorem 4. The solution of the problem i{36\)-08\) satisfies the following a priori 
estimate in Coo^)- 

\\u{x,t)\\^<\\u{x,0)U+ [ \\f{x,e)\\^de, (65) 

Jo 

whereas the solution of the problem 07l)-d39l) satisfies the a priori estimate 

\\u{x,t)\\,< ||n(a;,0)||i+ f \\f{x,e)\\^de. (66) 

Jo 

The estimates (1651 ). (I66l) complement the above a priori estimates (|60l)-(l62l) in 
the Hilbert spaces C2{^) and W^{n). 

3. Discretization in space 

In numerical solving transient problems, firstly we construct discretization in 
space. The resulting operator-differential equation should inherit the basic proper- 
ties of the differential problem, i.e, we speak of the positiveness (non-negativity) 
and self-adjointness of the diffusive transport operator as well as the adjointness of 
the convective transport operators in the corresponding finite-dimensional Hilbert 
spaces. We consider the standard finite difference approximations for model un- 
steady convection-diffusion problems in a rectangular domain. In addition, we 
discuss the problem of constructing approximations by means of the finite vol- 
ume method. 

3.1. Difference operators 

In a rectangle fi, we introduce a uniform in each direction grid. For grids in 
individual directions a;^, a = 1, 2, we use notation 
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where u stands for the set of interior nodes. On the set of grid functions that are 
equal to zero on the set of boundary nodes du (w = x ^2 = ujyjdoj), we define 
the Hilbert space H = ^2(0;) with the following scalar product and norm: 



{y,w) = ^y{x)w{x)hih2, \\y\\ = {y,yf'^ 



We use the standard index-free notations of the theory of difference schemes ll27ll . 
For the backward difference derivative, we have the representation 

U^= T • 

h 

Similarly, for the forward difference derivative, we get 

Ui+i - Ui 
Ux = -J • 

h 

Using the three-point stencil (nodes Xj, Xj+i), we can employ the central 
difference derivative 

Ui+i - 

Uo = ; . 

2h 

For the operator of the second derivative, we have 

Ux - Ux Ui_i - 2Ui + Ui+i 



The 2D difference diffusive transport operator is represented as the sum of the 
ID ones: 

2 

D = Y,D^^\ = a = 1,2, xGu;. (67) 

For smooth diffusion coefficients, we can assume 

a^^\x) = k{xi — 0.5/ii,X2), a^'^\x) = k{xi,X2 — 0.5/i2)- 



Properties of the elliptic difference operator D are well-known [|27L l25|l . For 
the 2D self-adjoint operator D, we have the lower bound 

D^D->±..E, M„ = | + |. (68) 
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We present also the upper bound for the diffusive transport operator, i.e., 

D < M3E (69) 

with the constant 

^3 = To 7; 

4 a^^\x) + a^^\xi,X2 + h2) 
+ 12 "^^^ n ■ 

Now we consider the difference analogs of the differential convective transport 
operators written in various forms. For the operator in the nondivergent form Ci, 
we put into the correspondence the 2D the difference convective transport operator 

2 

Ci = E^f"^' C'i"^Z/ = &^V", «=1,2, xeuj. (70) 

a=l 

In the simplest case of smooth enough convective transport coefficients, we as- 
sume 

b^°'\x,t) = Va{x,t), Xeu. 

Similarly, for the approximation of the convective transport operator in the 
divergent form C2, we employ the difference operator 

2 

C2 = 5]ci"\ C^'^)y=(6(")y)o^, a = 1,2, xeu. (71) 

a=l 

The approximation of the 2D convective transport operator in the skew-symmetric 
form is based on the representation Co = 0.5(Ci + C2) such that 

a=l ^ ^ 

a = 1, 2, X & Lo. 

Lemma 1. The difference operators Ca, a = 0, 1, 2 have the following properties 
of adjointness: 

d = -C2, C* = -Co (73) 
in the space of grid functions H. 
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Proof. It is easy to check directly that the ID convective transport operators in the 
divergent and nondivergent forms are adjoint to each other up to the sign. Taking 
this into account, we have 

= -{y,C2w) = iy,Clw). 

The skew-symmetric property of the operator Co follows from its definition as the 
half- sum of the operators Ci and C2. □ 
Similar properties can be proved for the 2D convective transport operators 
constructed with the use of coefficients Va{x,t), a = 1,2 shifted by a half-step 
in the corresponding directions. Such staggered grids are in common use in com- 
putational fluid dynamics 11351 l22n. 




Figure 1 : Control volume: □ — node for vi {x),t; O — node for V2 {x, t) 



Let us refer the convective-transport coefficient vi{x,t) with respect to the 
variable xi to the nodes of the grid which is shifted by a half-step along this 
direction. The grid for the coefficient V2{x,t) is shifted along X2 by 0.5/i2 (see 
Fig.ffl). 

For the difference convective transport operator in the nondivergent form, we 

get 

C^i'V = \ih'^'^ {xi - 0.5/11, X2, t)/^ + (xi + 0.5/11, X2, t)y^'). 
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2 

Ci = ^ci"\ xeu. (74) 

a=l 

For the difference convective transport operator in the divergent form, we em- 
ploy the representation 

C2y = Ciy+(b'^l^ +b^^Ay, xeu. (75) 

The following notation is used here for the difference derivative of the grid 
function given at half-integer nodes: 

_b{x + 0.5h) - b{x - 0.5/i) 

0« = ; . 

h 

This expression is a difference analog of the differential equality 

C2U = Ciu + div V u 

with a special approximation of div v. 

For the skew- symmetric convective transport operator Co = 0.5 (Ci + C2), 
from (l74l), (175]), we obtain 



1 



^0 1/ = TTT-^^'Ha;! + 0.5h,X2,t)y{xi + hi, 
Zri\ 



X2) 



- TTT-^^^Ha;! - Q.hhi,X2,t)y{xi - hi,X2), 
/All 

*^0^^^ = Wr^^^^^^^'^'^ + 0.5/l2,t)l/(Xi,X2 + h2) 
1112 

- 77^&''^Ha;i,X2 - 0.5h2,t)y{xi,X2 - /i2), 

Ztl2 

2 

a=l 

For the convective transport operators Ca, a = 0, 1,2, defined by (|74l)-(|76l). 
Lemma [H holds. 

In the multidimensional case, the inequality 

\{C^y,y)\<Mi\\yf, a = 1,2 (77) 
23 



with a constant Mi, independent of the grid steps, is also valid for the convective 
transport operators under consideration. For operators (TTOl) . dTTI ). we obtain a 
constant Mi, which depends on the first derivatives of the convective transport 
coefficients, whereas in the case (l74l) . (1751) it depends on the divergence, as in the 
continuous case. We formulate the corresponding statement using the following 
notation for the grids in separate directions: 



1,2, AT^-l, Naha = L}, 



1,2, ...,iV„, 



a = 1,2. 

Lemma 2. For the difference convective transport operators C^, a = 1, 2„ de- 
fined according to ( f7^ . ( fZTT ). f/ze estimates ([771) ^c/J w«Y/z constant 



Ml = - max max I fe^V I H — max max 16^^ , , 

2 aJiGwi ^2e<^f 



^(2)| 



2 xiea;+ 2:261^2 

whereas for f f7?l) . f f75l) — w/Y/i f/ie constant 



(78) 



Ml = - max 

2 tcGo; 



+ 

X2 



(79) 



Proof. Taking into account (1741) . (1751) . we get 



1 



{Ciy.y) = ~{C,y,y) = -{{Ciy^y) + {y^Cly)) 



(1) 



X2 



This implies the estimate (1771) . (1791) . 

For the operators (TTOl ). (1711) . we use the corresponding estimates for the ID 
operators, too. For instance, for the 2D convective transport operator in the non- 
divergent form, we have 

|(Ciy,y)| = |((Ci«+cf))y,y) 



< ^ max max \b'i^\{y,y) + ^ max max \b''^^\{y,y) 



(2), 



2 xi€uj+ X2&UJ2 



2 xiSi^i a-jgi, 



i. e., we arrive at the desired estimate (1771) . (fTSl) . 



□ 
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Finally, let us consider the subordination property of the convective transport 
operator to the diffusive transport operator under the standard restrictions k(x) > 

Ki > 0. 

Lemma 3. For the 2D convective transport dijference operators C = Ca, a = 
0, 1, 2, the following estimates hold: 

\\Cyf <M2{Dy,y), (80) 

where the constant M2 for the operators Ci, defined according to dTOl) . i^74\l . is, 
respectively, 

2 

M2 = —maxmax{b'-"\x,t)y, 

2 r 

M2 = — max < max (^'■^■'(xi — 0.5/ii, X2, t))^, 

max {h^^\xi,X2-^.hh2,t)f 



for the operators ([77]), ([^ 



M2 = — (2 max max {h'^'^\x,t)f + Momax (h^P + hP\] 

Ki \ a xeuj ^ ' xew \ xi X2 / J 



M2 = — 2 max < max (6(')(a;i -0.5/11, X2,t))' 

1^1 \ I xGui^ XUI2 



max (6(2)(a;i,x2-0.5/i2,t))' 



+Momax(bP +b?A ] , 

x€u} \ xi X2 / J 

and for the operators ^72\l . ^76ti — 

M2 = — (3 max max (b'^''\x,t)Y + Momax (b^ + b^^ ^ 

Kl y a xeu! ^ x^u) \ xi X2 

M2 = — 3 max < max (6(')(a;i -0.5/ii,X2,t))', 
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max 

XEUJI XUJ2 

+ Mo max + ) , 



where Mq — is the constant from ( A68\l . 
Proof. Taking into account the inequality 



9 \ 2 p 



for operator ( TTOl ). we have 

XdU) 

< 2maxmax(6(°)(a;,t))2||V?/f < M2{Dy,y). 



The estimate for the operator (|74l) is obtained in the same manner. 
For the discrete operator dTTI) . we employ the representation 

C2y = ^ {b^'\x, - h,,X2,t)y^' + h^'\x^ + h,,X2,t)y^' 
+ b^^\x^, X2 - h2, t)y^' + b^^^ (xi, X2 + h2, t)y''') 
+ fen 6^5) 

\ Xl X2 / 



Thus, 



\\C2yr < 2 5^(6«(xi - h,X2)fiy'^Yhh2 
+ 2j2{b^'Hxi + hi,X2)fiy'^fhh2 

X£U! 

+ 2j2ib^'\^i,X2-h2)fiy''fhh2 

XdU) 
Xduj 

+ 2max(6i') + 6i'))2(y,y). 

CKGOJ X\ X2 
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In view of the Friedrichs inequality, we obtain the estimate (f80l ) with the constant 
M2 given in the lemma. 

For the grid operator in the divergent form with coefficients on staggered grids, 
on the basis of (l75l) . we have 

2 



|C2yr = 2||Ciyr + 2 



y 



For the first term in the right-hand side, we use the already derived estimate for 
Ci, whereas for the second one we apply the Friedrichs inequality. 

Subordination estimates for the difference 2D convective transport operators in 
the skew- symmetric form Cq are derived from the estimates for the operatorsCo, a 
1,2. □ 

The above values for the subordination constant M2, in spite of their awkward- 
ness, demonstrate the fundamental independence of this constant of the compu- 
tational grid. The constant M2 depends on the values of the convective transport 
coefficients Va{x,t), a = 1,2 (the velocity) and on div v, to be more precise, on 
their difference approximations. 

In numerical solving the problem (|4TI) . (|59l ), using the above discretization in 
space, we obtain the operator-differential equation 

^ + Ay = ip{x,t), A = C + D, xeuj, 0<t<T, (81) 

y{x,0) = u\x), xeu. (82) 

To investigate this semi-discrete problem, we employ the above properties of the 
difference operators of convective and diffusive transport. In particular, we can 
obtain analogues of a priori estimates of Theorem [2l 

3.2. Monotone schemes for 2D convection-diffusion problems 

For 2D difference convection-diffusion problems in nondivergent and diver- 
gent forms, the maximum principle is formulated. Unconditionally monotone 
schemes are constructed on the basis of the regularization principle for difference 
schemes. 

To simplify the material presentation, we will consider difference schemes for 
stationary 2D convection-diffusion equation on uniform rectangular grids. The 
corresponding discrete analogues on the standard five-point stencil cross are writ- 
ten in the following form: 

'y{x)y{x) - ai{x)y{xi - /ii, X2) - Pi{x)y{xi + hi, X2) 

- a2{x)y{xi, X2 - h2) - l32{x)y{xi, X2 + h2) = ^{x), x e u. 

(83) 
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This grid equations are considered with the boundary conditions 

y{x) = 0, xe duo. (84) 
Assume that the coefficients of the difference scheme (|83l) satisfy the conditions 

> 0, > 0, j = l,2, 7(0;) > 0, xeu (85) 

and suppose that 

aj(a;) > 0, /3j(a;) > 0, j = l,2, x G 

We highlight two classes of the monotone difference schemes ([83] ). (f84l ). i.e., 
the schemes that satisfy the difference maximum principle. 

Theorem 5. Assume that in the scheme d^?l)-([S5l) ^{x) > for all x ^ u (or 
(p{x) < Oforx G uj). Then for 

'^(x) > ai{x) + a2ix) + (3i{x) + I32{x), xeu (86) 

or for 

l{x) > ai{xi + hi,X2) + - hi,X2) ^^^^ 
+ a2{xi,X2 + 112) + l32{xi,X2 - 112), X e 00 
we have that y{x) > Ofor all x E u (y{x) < Ofor all x G u). 

Proof As usual, the argument is by reductio ad absurdum. Suppose, e.g., that the 
conditions (f86l ) are valid and the solution of equation ([83] ) with the non-negative 
right-hand side is not non-negative at all grid points. Let x* be an interior grid 
point, where the solution has the minimal negative value. If this value is achieved 
at several points, then we choose the grid point, where y{xl — hi, x^) > y{x*). 
We write equation (|83l) at this point as 

'j{x*)y{x*) - ai{x*)y{xl - hi,xl) - (3i{x*)y{xl + /ii, x^) 
- a2{x*)y{xl, x*2 - 112) - (32{x*)y{xl, x*2 + h2) 
= (p{x*), X* G u. 

Under the theorem conditions, the right-hand side is non-negative, whereas the 
left-hand side, in view of ([85]), (l86l) . 

(7(a;*) - ai{x*) - f3i{x*) - a2{x*) - f32{x*))y{x*) 

+ aiix*){y{x*) - y{xl - hi, x^)) 

+ (3i{x*){yix*)-y{x*i + hi,x*2)) 

+ a2{x*){y{x*) - y{xl, x^ - /la)) 

+ P2ix*)iyix*)-y{xl,x*2 + h2))>0. 
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We arrive at a contradiction, and therefore y{x) > for all a; G cu. 

Now we consider a more complicated case of the difference scheme (|83l) . (f84l) 
with the conditions ((l85l). (l87l) satisfied. Suppose that for the non-negative right- 
hand side of equation (l83l) . there exists a subset of the interior grid points to*, 
where 

y{x) < 0, a:; G w*. 
Summarize equations (l83l) over all these nodes: 

(^7(a;)y(a;) - «i(a;)?/(xi - hi, X2) - (3i{x)y{xi + hi, X2) 
- a2{x)y{xi,X2 - h2) - f32{x)y{xi, X2 + /i2)) 

ccgo;* 

For the left-hand side of this equality, we have 

(lix) - ai{xi + hi, X2) - - hi,X2) 

- a2{xi,X2 + h2) - (32{xi,X2 - h2)^y{x) 

+ X] ("1(^1 + ^2)1/(3^) - ai{x)y{xi - hi, X2)) 

+ ^ i(3i{xi - hi, X2)y{x) - Pi{x)y{xi + hi, X2)) 

Xdu)* 

+ X] («2(a;i, + h2)y{x) - a2{x)y{xi,X2 - /i2)) 

X(^UJ* 

+ ^ (/32(a;i, - h2)y{x) - /32(a;)2/(a;i, 0:2 + /i2))- 

We see immediately that each of these terms is non-negative and they cannot 
be equal to zero simultaneously. Thus, we again obtain a contradiction. □ 

The maximum principle for multidimensional difference equations, where the 
sufficient conditions of type (f86l) are satisfied, is well-known in the theory of dif- 
ference schemes [27]. For elliptic difference problems, the maximum principle 
in the form (l87l) is presented in work [IjtIi . Similarly to the ID case, the condi- 
tions (|86l) may be associated with the condition of diagonal dominance by rows 
if we use the natural order for the points of the numerical solution, whereas the 
conditions ([87] ) correspond to diagonal dominance by columns. 
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Consider the 2D boundary value problem for the convection-diffusion equa- 
tion in nondivergent form: 

a=l a=l ^ ' 

u{x) = 0, dn. (89) 

For numerical solving the problem (l88l) . (|89| ). we employ the scheme 

Ciy + Dy = <^{x), xeu. (90) 



We restrict ourselves (see [|32n for details) to the 2D operator of convective trans- 
port in the form 

V = ^rr\ik{^i - 0.5/11, x,)y^^ + k{x, + 0.5/11, a;2)y"^), 
Cf\ = ^^{k{x^,X2 - 0.5/i2)/^ + k{x^,X2 + 0.5/i2)z/"^), 

2 

C^ = J2c["\ xeu. (91) 

a=l 

We formulate the conditions of monotonicity for the difference scheme (l90l) . 
(|9T]). We write it in the form ([83]), (f84l) with 



/ Vi(x) 1 \ , / 



7(a;) = ai{x) + a2ix) + f3i{x) + /32{x), xeu. 
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The monotonicity condition (f87l) is obviously satisfied, and the positiveness of the 
coefficients of the difference scheme ([83] ) (the condition ([85] )) leads to the natural 
restrictions 

\daix)\ < I, a = 1,2, xEu, 

where 

Unconditionally monotone difference schemes for the 2D convection-diffusion 
equation ([88]) . ([89]) are constructed by means of regularization (disturbance) of the 
diffusion coefficient. Instead of ([90] ). we consider the difference scheme 

2 

Ci|/ + 5^(l + f?,)D(°)i/ = V^(a;), xeu. (92) 
The scheme ([92]) is monotone under the condition 

1 + Qaix) > \9a{x)\, XEU, a = 1,2. 

We present some variants of regularizing grid functions Qaix), a = 1,2,, which 
lead to unconditionally monotone difference schemes ([92|) . 
For example, the choice 

1 + Qaix) = 6a{x)cth.6a{x), X ^ U, = 1,2 

corresponds to the use of exponential schemes [0, Isl] for each individual direc- 
tion. It is also possible to select the scheme with 

Qaix) = r]9l{x), X eu, a = 1, 2, 

which are monotone if i] > 0.25. Among unconditionally monotone difference 
schemes, we highlight the regularized scheme ([92l) . where 



6'^(x) 

^"(^) = 1 , T/3 / M ' X eu, a = 1,2. 

1 + \9a{x)\ 

The scheme with the upwind differences corresponds to 

Qaix) = r]\6a{x)\, X E UJ , = 1,2 
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with ?7 = 1. In this case, the difference operator of convective transport seems like 
this: 

(1) k{xi-0.5hi,X2) +, . , k{xi + 0.5hi,X2) , 
k{x) k[x) 

(1) A;(xi,X2 - 0.5/12) . ^2 , A;(xi,X2 + 0.5/12) . 

2 

C^ = Y,C["\ xeu. (93) 

0=1 

Obviously, the corresponding scheme has only the first-order approximation. 
Now consider the convection-diffusion equation in the divergent form: 



x:^m-)u)-x:^(m-)|^ 

a=l a=l ^ 



fix), xen, (94) 



which is supplemented with the boundary conditions (|89l) . For this equation, we 
consider the difference scheme 

C2y + Dy = (p{x), xeu. (95) 

Define the difference operator of the convective transport as 

= T^M^i + 0.5hi,X2){y{xi + hi,X2) +y{x)) 
Zri\ 

~ 7TT-^i(a;i - 0.5/ii,X2)(y(xi - hi,X2) +y{x)), 
2hi 

^2 y = 7T7-^2(a;i,a^2 + 0.5/;,2)(?/(a;i,X2 + /i2)+y(a;)) 
2h2 

- 777-^2 - 0.5/2,2) (^(xi, X2 - h2) + y{x)), 
2hi 

2 

C2 = ^C("\ xeu. (96) 

a=l 

This is typical for coefficients of convective transport defined on staggered grids. 
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The difference scheme (|95] ). (|96l ) may be written in the form ([83] ) with the 
coefficients 

tiifxi — 0.5/ii, X2) A;(xi — 0.5/ii, X2) 
"i(^) = + n ' 



2hi hi 



2/ii /if 



, , t;2(xi,X2 - 0.5/12) /cfxi, X2 - 0.5/12) 
= 2h, + hi ^ 

^2(3^1,3^2 + 0.5/12) , /c(xi, X2 + 0.5/12) 



= 2)1^ + hi 

^{x) = ai{xi + hi,X2) + Pi{xi - hi,X2) 

+ a2{xi,X2 + h2) + l32{xi,X2 - h2), xeu. 
Therefore, the condition ( [87] ) is valid, and from ([85] ). we get 

|6'i(a;i — O.5/11, ^2)1 < 1, X E Ui X U2, 

\92{Xi,X2 — 0.5/12)1 < 1, X E Ui X UJ2 ■ 

A class of regularized difference schemes for the convection-diffusion equa- 
tion in the divergent form is defined as 

C2y - ((1 + gi{xi - 0.5hi,x2))k{xi - 0.5hi,x2)y''^)x^ 

- ((1 + f?2(xi, X2 - 0.5h2))k{xuX2 - 0.5/12)/^)., (97) 
= ip{x), X E u. 



Sufficient conditions for the monotonicity of the scheme ([97] ) are written as 

1 + Qi{Xi — 0.5/11,^2) > \Oi{Xi — 0.5/li,X2)|, X E Ui X U2, 

1 + ^2(a;i,a;2 - 0.5/12) > 1^2(3^1,3^2 - 0.5/i2)|, x E uJi x w^. 

Some approaches to select regularizing grid functions Qa{x), a = 1,2 in 
order to obtain unconditionally monotone difference schemes were considered 
above for the convection-diffusion equations in the nondivergent form. 
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3.3. Triangular grids 

The possibilities of solving boundary value problems for PDEs on irregular 
grids are discussed here. The focus is on constructing difference schemes on 
triangular grids (as the most common unstructured grids). We emphasize approx- 
imations on Delaunay grids (triangulations) that demonstrate optimal properties. 
The problem of discretization in space is illustrated considering steady- state prob- 
lems. 

The basis for the construction of discrete analogs is the balance method (the 



integro-interpolation method) [I26IJ27!]. which nowadays (in the English literature) 
is referred to as the finite volume method [42, 15]. The efficiency of this approach 
becomes evident in designing difference schemes on irregular grids. As a control 
volume in Delaunay triangulations appears Voronoi diagrams. 

Among general irregular grids we distinguish structured grids that are topo- 
logically equivalent to regular grids. A typical example of unstructured meshes 
are triangular grids. There are discussed general issues of designing grids and 
discretization on them. 

Numerical solving boundary value problems of mathematical physics in com- 
plicated domains is carried out using irregular grids. A computational domain Q 
is assumed to be irregular (nonrectangular and not composed of rectangles). Be- 
cause of this, we have to use nonrectangular grids. Among irregular grids, we 
emphasize two main classes. 

Structured grids. The most important example of such a type of grids is irreg- 
ular quadrangular grids that inherit, in many senses, properties of standard 
rectangular grids (they are topologically equivalent to them). 

Unstructured grids. In this case, a stencil of a difference scheme has no fixed 
structure. It is impossible to connect topologically such a computational 
grid to a regular rectangular grid. In particular, schemes have a different 
number of neighbors at each grid point. 

Approximations on structured grids can be performed on the basis of the 
above-mentioned closeness of these grids to standard rectangular grids. The sim- 
plest realization of this situation is to use new independent variables. In this case, 
a grid that is irregular in the original coordinates is transformed into a regular one 
in new independent coordinates. 

The second possibility is not associated with a formal introduction of new 
coordinates; it is implemented using an approximation of the original problem on 
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an irregular grid. It is clear that the use of simple approaches for the construction 
of difference schemes on the basis of uncertain coefficients for irregular grids is 
possible, but it seems not so reasonable. 

Advantages of structured grids results from the conservation of the canonical 
structure of neighbors for each grid node, i.e., the conservation of the stencil. 
This simplifies, in particular, the process of programming and solving difference 
problems. But the problems of constructing difference schemes on such grids are 
not much less difficult than for the general unstructured grids. 

Among structured grids, it is necessary to distinguish an important class of 
orthogonal grids. In this case, the advantages of structured grids over unstructured 
ones become evident because a lot of problems connected with the development of 
difference schemes, and the solution of grid equations is radically simplified. If it 
is necessary to use the advantages of structured irregular grids over unstructured 
ones, it is better to restrict ourself to orthogonal curvilinear grids. Problems of 
grid generation are not necessarily more difficult than the problems being solved. 
Moreover, this situation is the most typical one. Therefore, it is better to make 
efforts (that are comparable to the solution of the original problem) to optimize 
the computational grid. In complicated computational domains, it is reasonable to 
use the multiblock technology of generation of orthogonal grids that is based on 
modern CAD systems. 

An arbitrary grid is generated from a set of nodes. The most simple and natural 
approach is to define a triangulation, i.e., to construct a triangular grid. There is 
no need to use more complicated structures of unstructured grids. 

For the given points, a triangulation can be performed in different ways. Note 
also that for a given set of nodes, we obtain the same number of triangles by any 
triangulation method. 

Thus, we need to optimize the triangulation by some criteria. The main opti- 
mization criterion consists in the following: the obtained triangles should be close 
to equilateral ones (they should be without too sharp angles). This is a local cri- 
terion governing to an individual triangle. The second (global) criterion declares 
that adjacent triangles must not differ too widely in an area — the criterion of grid 
uniformity. 

There is a special triangulation — the Delaunay triangulation [[ill 0], which 
has a number of optimal properties. One of them is the tendency of obtained tri- 
angles to be equiangular ones. The above-mentioned property can be formulated 
more exactly in the following way: in the Delaunay triangulation, the minimal 
value of inner angles of triangles is maximized. The formal definition of the 
Delaunay triangulation is associated with the property that for each triangle all 
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the other nodes are located outside of the circumcircle. For our further presen- 
tation, the relation between the Delaunay triangulation and the Voronoi diagram 
(tesselation) is very important. 

A Voronoi polygon for a separate node is a set of points lying closer to this 
node than to all the other nodes. For two points, the sets are defined by the half- 
plane bounded by a perpendicular to the middle of the segment connecting these 
two points. The Voronoy polygon thereby will be the intersection of such half- 
planes for all pairs of nodes created by this node and all the other nodes. Note that 
this polygon is always convex. 

Each vertex of a Voronoi polygon is a point of contact of three Voronoi poly- 
gons. The triangle constructed by the corresponding nodes of contacting Voronoi 
polygons is associated with each of these vertices. This is exactly the Delaunay 
triangulation. Thus, between the Voronoi diagram and the Delaunay triangulation 
a unique correspondence is established. 

In the case of the Delaunay triangulation, we obtain the optimal decomposition 
of a computational domain according to the given set of nodes. The decomposi- 
tion is optimal in terms of maximization of minimal angles of triangles. For the 
Delaunay triangulation, there does exist the corresponding Voronoi diagram that 
uniquely determines a set of points of the domain for each node. This separation 
of the set points is made by the clear geometrical criterion of optimal closeness 
to the node. Thus, the Delaunay triangulation and the Voronoi diagram determine 
completely (optimally and uniquely) a computational triangular grid and a control 
volume. 

The Delaunay triangulation is widely used in numerical practice for construct- 
ing finite element procedures. There also exist a lot of developed numerical meth- 
ods for generating such triangular grids, the appropriate software is also available. 

3.4. Dijference schemes on triangular grids 

We start with a discussion of some possible general approaches that may be 
applied to (and find) practical applications. 

The simplest (from the methodological viewpoint) approach to construct dis- 
crete analogs on triangular grids consists in using the finite element procedures 



[|34ll . However, it is not always possible to employ standard variants of finite ele- 
ments. 

In computational practice, there are widely used piecewise linear finite ele- 
ments, which correspond to the approximation of the numerical solution on each 
triangle via linear functions. In the convection-diffusion problems, we obtain ana- 
logues of schemes with central difference approximation of the convective terms. 
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In constructing finite element approximations, there do exist some problems to 
obtain monotone procedures, i.e., the schemes that satisfy the maximum principle. 
In the theory of finite elements, the problem is resolved not using the Galerkin 
method, but using its generalization — the projection Petrov-Galerkin method. In 
this case, probe functions that are used to construct the solution, and test functions 
that generate a system of equations, are distinct. In this sufficient artificial way, 
we obtain schemes, which are, e.g., very similar to the schemes with the upwind 
differences 

In the method of support operators fl^, the original problem is formulated in 
terms of differential operators from the vector analysis: the divergence, gradient, 
and rotor. Next, only one of them is freely approximated on a selected grid. Other 
operators are defined by some prescribed relations of integral nature that exist 
between the differential operators. This ensures consistent approximations of the 
operators that provide the fulfillment of such essential properties as conservatism, 
adjointness and so on. 

The method of support operators has been developed by many researchers 
just for triangular grids. The main peculiarities are associated with a selection 
of the set of grid functions. For example, the solution can be approximated at 
vortices of the triangular grid, whereas fluxes can be refereed to cells (the cell- 
vortex arrangement of variables) or to the the midpoints of cell faces (staggered 
grids). 

One of the basic approaches to the construction of difference schemes on irreg- 
ular grids is the classical integro-interpolation method [27]. This balance method 
is based on the following main points. First, we must specify a grid (determine 
a set of nodes and a set of grid functions). Secondly, for each node, we define a 
neighbor domain (control volume) — a part of the computational domain adjoin- 
ing to a given computational node. And finally, a difference scheme is obtained by 
integrating the original equation over the control volume using some assumptions 
about the solution behavior. A set of these three components specify a particular 
variant of the control volume method. 

In constructing difference schemes on triangular grid via the control volume 
method, it seems natural to set grid functions at the grid nodes. This is the stan- 
dard, but not the only variant. As an alternative, we can reffer grid functions 
to some points being connected with a triangular cell. The specification of grid 
functions at the vertices of the Voronoi polygons provide an example. 

The second problem is connected with selecting a control volume. During a 
triangulation procedure, in many approaches, a part of the triangle appearing from 
the intersection of medians is separated as as a control volume. In this case, each 
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node obtains a part of the triangle with an equal area. 

An interesting variant of the control volume selection is associated with the 
Voronoi diagrams. In this case, each individual node has a part of the whole 
computational domain that is closest to it. In this case, there is no division into 
triangles with equal areas. It is important that in the both above approaches, a 
criterion for selecting a control volume is clear and is connected with geomet- 
ric requirements: in the first case — the triangle is decomposed into three parts 
of equal area, in the second — we obtain geometric proximity of points of the 
computational domain. 

Among the merits of the division by medians, we emphasize that it may be 
conducted for an arbitrary partitioning into triangles, i.e., not only for the De- 
launay triangulation. Advantages of the Voronoi tesselation are more essential 
and associated with the orthogonality of the triangle sides to the faces of Voronoi 
polygons. 

Heuristic arguments in favor of Voronoi polygons are associated with the idea 
of globalization (optimization) of grids and control volumes — optimization for 
all nodes, rather than for a single triangle. 

The balance method for triangular grids need to be implemented using the De- 
launay triangulation with Voronoi polygons as control volumes. This is the most 
natural way that allows to construct difference schemes with optimized triangular 
grids H- 

3.5. Diffusive transport operator 

Assume that a computational domain is a convex polygon Q with the boundary 
dil. The points of the domain are denoted hy x = (x^^\ x^"^^). 

In the domain 1] = i7 |J d^l, we consider the grid ZJ, which consists of nodes 
Xi, i = 1,2, . . . , M, and the angles of the polyhedron are nodes. Let w be a 
set of interior nodes and dcu is a set of boundary nodes, i.e., cu = iUf]il, du = 

mf]dn. 

Each node Xi, i = 1,2, . . . , M, is associated with a certain part of the com- 
putational domain i7j treated as a control volume. A Voronoi polygon or its part 
belonging to are selected as the control volumes. A Voronoi polygon for an in- 
dividual node is a set of points lying closer to this node than to all the other ones. 
For two nodes, the sets are defined by the half -plane bounded by the perpendic- 
ular to the midpoint of the segment connecting these two nodes. The Voronoi 
polygons thereby will be the intersection of such half -planes for all pairs of nodes 
created by this node and all the other nodes. Each vertex of a Voronoi polygon 
is a point of contact of three Voronoi polygons. The triangle constructed by the 
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corresponding nodes of contacting Voronoi polygons is associated with each of 
these vertices. This is exactly the Delaunay triangulation. 

Control volumes cover the whole computational domain, so that 



M 



1=1 



i^j, i, J = 1,2, M. 
For the common faces of control volumes, we use notation 

dnif]dnj = Tij, i^j, z, j = 1,2, M. 

For the node i, we define a set of neighboring nodes >V(i) that have the control 
volumes with common faces with the control volume for the node i, i.e., 

>v(i) = {j I an, fi on, j = 1,2, m}, 1 = 1,2, ...,m. 

Litroduce notation 



V, 



cii r 



dx, i, j = 1,2, M, 



for the area of the control volume and the length of the edge of the Voronoi poly- 
hedron, respectively. 

For the grid functions y{x), w{x) that are specified at the nodes x & u and 
vanish at the boundary nodes x e du, inH = L2{u), we define the scalar product 
and norm 

iy^w) = ^ Viy{xi)w{xi), \\y\\ = {y,vf''^. 
Define the distance between the nodes 



a=l 



1/2 



and the midpoint of the segment connecting these nodes as follows 



X 



^3 



(a) 



(a) . (a)\ -I r, 

XI ^ + x) a= 1, 2. 
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For simplicity, we assume that the coefficients of the convection-diffusion 
equation and its solution itself are sufficiently smooth. The discrete operator of 
diffusive transport corresponding to the interior node of the computational grid 
Xi E cu is defined according to the integro-interpolation method by means of the 
integration over the control volume fi^: 



Du = {Du)i ~ ^ / Vudx. 



(98) 



For the diffusive flux vector, we use the expression 

q = —k{x) gradw, 

so that 

Vu = div q. 

From this, we obtain for the right-hand side of (|98]) that 

J Vudx = j{q,n)dx, (99) 

where n is the outer normal. 

A difference approximation for the normal component of the diffusive flow 
through the face 'jij is denoted by q'lj, and therefore, from (|98] ). (|99| ). we get for 
the difference diffusive transport operator the following representation: 

' ie>v(i) 

In the case of smooth enough coefficients and the solution itself, it is natural 
to employ elementary approximations for the flux along the normal at the point 

Xij : 

«S = -*(-«)|;=|)- 

From (1 1001) . (IIOII) . we obtain the difference operator of the diffusive transport: 

{Dy)i = -^ k,k{Xij) ^'~y XiEu. (102) 
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As in the case of difference schemes on regular grids, it is possible to introduce 
various approximations for the diffusive flows. This issue is very important, in 
particular, for problems with piecewise-smooth coefficients of the equation. 

For the grid functions yi = y{xi) = 0, Wi = 0, Xi E dco, we have 

This summation over all faces of the Voronoi polygons for all interior nodes may 
be rewritten in the more convenient form: 



^ ^ ^ ^ lijd^Xi, Xj)k{xij] 



yj yi Zj Zi 



Thus, {Dy, w) = {y, D*w), i.e, the difference operator (11021) is self-adjoint in H. 
In view of 

iDy,y) = lYl E ' ^^^^^ 

it is also positive (D = D* > 0). 

Now we establish a discrete analog of the Friedrichs lemma. In various for- 
mulations, it is proved in the theory of finite elements for considering difference 
schemes on irregular structured and unstructured grids. Difference schemes that 
are based on the Delaunay triangulation and the Voronoi diagram demonstrate 
some peculiarities, and therefore the discrete analog of the Friedrichs lemma must 
be proved for them separately. 

Lemma 4. For the grid functions y^ = y{xi) = 0, Xi G du, the following in- 
equality is valid: 



with the constant 

Mo = 1 + 1, 
° 16 16' 

where la, « = 1,2 are the side lengths of the rectangle with sides parallel to the 
coordinate axes that contains completely the whole polygon Q. 
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Proof. On the set of the grid functions yi = y{xi) = 0, ajj G duo, in view of (11021) . 
we define the discrete Laplace operator A via the expression 

1 , -y' 



Using this notation, the inequality (|104l) may be rewritten in the more compact 
form: 

||yf <Mo(Ay,y). (106) 

To estimate the lower bound of the discrete Laplace operator, we employ the 
solution of an auxiliary boundary value problem. We will show that in the in- 
equality (11061 ) under consideration, we can put 

Mo = maxty(a;), (107) 

where w{x) is the solution of the problem 

Aw = 1, xeu. (108) 

We consider in H the eigenvalue problem 

Ay = Xy, xeu (109) 

for the grid operator A = A* > 0. For the problem (11091) , the following inequality 
holds: 

{Ay,y)>Knin\\yf (HO) 

for any y{x). The equality in (II 101 ) is achieved only for the eigenfunctions v(x) 
that correspond to the minimal eigenvalue Amin- Thus, the estimate (11061) will be 
established if we will show that Mq^ < X^-m. 

First of all, let us explain that v{x) is a constant-sign function. Suppose that 
this is not true and v{x) changes its sign on the grid to. Now consider the function 
|t>(a;)|. Then taking into account 

we obtain 

{A\v{x)\,\v{x)\) ^ {Av{x),v{x)) 
i\v{x)\, \v{x)\) {v{x),v{x)) ■ 
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This contradicts the fact that the minimal Rayleigh ratio 

is achieved for y{x) = v{x). 
From Av = XmmV, we have 

{v{x),i) ■ ^^^^^ 

In view of (11081) . we obtain for the denominator 

{v{x), 1) = {v{x), Aw{x)) = {w{x), Av{x)) < maxw{x) {Av{x), 1), 

since Av{x) > 0, x E u. Thus, from (lllll) . it follows immediately that M^^ < 
Aminj and for the constant Mq, we can use the expression (11071) . 

Next, we apply the maximum principle to discrete elliptic equations. We place 
the polygon into the rectangle 

no = {x = (a;(^), I a„ < x^^^ < b^, a = 1, 2}, 

where L = 6„ 1,2. Consider the function 

W{x) = -/i((a;« - ai)(a;(^) - h) + {x^^^ - a2){x^^^ - 62)) 

with some positive constant fi. We write this majorant function as 

W{x) = - x\'^ f + - V) + gi{x), 

where gi{x) is a linear function. 
For linear functions, we have 

Agix) = 0. 

To show this, it is sufficient to consider the function g(x) = x'^^\ By virtu of 
(I105D . we obtain 

(1) (1) 

1 , X) ' — X) 



Vi ^ ^ d(xi,x.i] 



We have 



xf^ - ^ 



43 



where (^ji is the angle between the segment [xi, Xj] and the axis Ox^^\ In this 
case, lij cos(v9ji) is the projection of the j-th face of the Voronoi polygon ^7j on 
the axis Ox^^\ For a closed polygon 



Aa;^^) = ^ lijCos{<^ji) = 0. 



Vi 

Taking this into account, for the function W{x), we obtain directly 

{AW{x))i = li{k(f{x, Xi))i = XI ^iA^i^ Xj) = 4- 

That is why for = 0.25, we have 

AW{x) = 1, X euj. 

At the boundary nodes 

W{x) > 0, X e du, 

and therefore the function W{x) is a majorant for the problem (|108|) . For the 
constant Mq in the inequality (11041) . the following estimate holds: 

Mo<maxW{x) = + z^. 
xeuj ib ib 

This completes the proof of the lemma. □ 
This makes possible to formulate the main result concerning the properties of 
the difference operator of diffusive transport. 



Theorem 6. For the grid operator of diffusive transport determined from m4\l . the 
inequality 

D = D*>—E (112) 
- Mo 

is valid for the grid functions from the space H = L2{uj). 

It is important that the constant Mq in the Friedrichs inequality (11041) is in- 
dependent of nodes of the computational domain, and the estimate itself is quite 
similar to the estimate for the differential operator. 
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3.6. Convective transport operators 

In the construction of difference operators of convective transport, we start 
with the operator in the divergent form (|44|) . Assume that 

C2U = C2UdX (113) 

for the interior nodes of the grid. For the right-hand side, we have 



C2udx = J {v, n)udx. 

Similarly to the case of the diffusive transport approximation, the normal com- 
ponent of the velocity is referred to the midpoint of the segment connecting grid 
nodes. Introducing notation 

bij = {v,n){xij), 
from (II 131) . the difference operator of convective transport is written as 

Now we discuss approximations of the convective transport operator in the 
nondivergent form (l43l) . There do exist several opportunities. The first one is con- 
nected with the use of a sufficiently complicated (not so evident) structure for the 
difference convective transport operator in the nondivergent form based on special 
formulas of approximate integration. But we have another way. Using the idea of 
the method of support operators, first, we design a simple approximation of the 
convective transport operator in the divergent form (II 141) . To obtain a difference 
approximation of the operator (|43l ). we search a difference operator that is adjoint 
to (|114l) . In doing so, we get a difference operator of the convective transport in 
the nondivergent form. Such an opportunity is essential for developing approxi- 
mations on irregular grids. 

Straightforward calculations yield 

{C2y,w) = ^J2 Yl ^iAjVi^i + ^iAjVi^i 

Xieuj jeW{i) Xieuj jeW{i) 

= -{y,Ciw), 
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where we take into account that 6^ = —bji. This allows to define the difference 
operator 

{Ciy)^ = — hjbij ^^ ^ ^\ XiEu. (115) 

* jew(i) 

By construction, the difference operators of convective transport in the divergent 
and nondivergent forms determined in accordance with (II 141) . (II 151 ) are adjoint to 
each other with within the sign, i.e., 

CI = -C2. (116) 

As for the operator of convective transport in the skew-symmetric form (|45l) . 
using the representation 

Co = l{C, + C2), 
from (11141) and (II 151) . we get the most compact approximation 

(Col/)i = ^ Y hjbijVj, XiEuj. (117) 

The primary feature of this difference operator is 

Co* = -Co, (118) 

and moreover, this skew-symmetric property is true for any velocity field — it is 
valid for arbitrary vectors v, not necessarily satisfying some difference analogue 
of the incompressibility constraint 

d-^-E^ = 0' (119) 

a=l 

To study difference analogs of the boundedness of the convective transport 
operator and its subordination to the diffusive transport operator, we discuss the 
important features of the approximations (II 141) and (II 151) in detail. 

For numerical solving continuum mechanics problems, it is essential to have 
consistent approximations of the convective transport operator in the divergent 
and nondivergent forms. The consistency is treated in the sense that one dif- 
ference operator coincides with other operator if the corresponding difference 
incompressibility constraint holds. This issue is very important due to the fact 
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that this equivalence takes place for the differential equations, and just it ensures 
the fulfillment of several conservation laws. In particular, for the incompress- 
ible Navier-Stokes equations, the convective transport operator in the momentum 
equation makes no contribution neither to the kinetic energy, nor to the individual 
momentum components (we speak of energy neutrality and neutrality with respect 
to the momentum). Unfortunately, elementary approximations ensure only one of 
these properties — either for the kinetic energy or for the momentum. 

For an incompressible fluid (the constraint (II 191) ). at the differential level, we 
can use any of the above-mentioned three equivalent forms of the convective terms 
(t43l) - (t45l) . It follows from the formula of vector analysis: 

dw{vu) = {v ■ gradu) + u divv, 

which, in turn, is based on the differentiation formula of the product of two func- 
tions. For the operators of convective transport, we have 

C2U = Ciu + diw u. (120) 

In constructing difference approximations for the convective transport opera- 
tor, it seems natural to develop such difference operators that satisfies the property 
(11201) of differential operators, i.e., the property of equivalence of various differ- 
ence approximations. 

From (11141) . (I115L we get 

C2y = Ciy + y ^ kjbijyi, Xi e u. 

Transform this equality to the form that is similar to (11201) : 

C2y = Ciy + divhV y, (121) 
where the difference operator of divergence is 

divhV = y ^ hjbij, XiEu. (122) 
* je>v(i) 

An approximation of the divergence for a vector by means of (11221) is obtained 
using the definition 

£fdS 

div V = lim 



5^0 J dV 

V 
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where V denotes the area, dV stands for the boundary of the domain, and S is 
the domain diameter. The expression (11221) may be treated as the correspond- 
ing quadrature formula for the right-hand side in the integration over the control 
volume for the node Xi E cu. 
In view of (11211) . we write 

Ciy = Coy-^ divh v y, 

{ (123) 
= Coy + - divh v y. 

In view of (II 181) . from (11231 ). it follows immediately that 

\{C^y,y)\<M4yf, a = 1,2, (124) 

with a constant Mi that depends only on the compressibility of a medium (at the 
discrete level), i.e.. 

Ml = ^||div^'«;||oo. (125) 



Here we use notation 



\y\\oo = max\y[Xi 



for the norm in Loo{uj). To prove (I124|) . multiply equation (11231) scalarly by y and 
apply the skew- symmetric property of the operator Cq. 

To obtain a difference analogue for the inequality representing the subordi- 
nation of the difference convective transport operator to the difference diffusive 
transport operator, we start with the upper bound for the expression 



2 



Xieuj \ieW(i) 



\\c^y\\' = J2v.\ E • (126) 

XieL 

For the right-hand side, we have 



2 

yj - yi 



Xieuj \jeyv{i) 



< max max — ( lijd(xi,xA^ — 



XiEui \jeyv{i) 
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Taking into account the inequality 



and assuming 



~ {hjd{xi, Xj)) ^ , — ihjd{xi, Xj)) ^ . 

dyXi^ Xj J 

we obtain 

2 



- X/ 4y: lijd{xi,Xj) ^ lijd{xi,Xj) (j ^ ^j ■ 
In view of 

je>v(i) 

substitution into (11261) yields 

\\Ciyf < max max V V l^A^^.^j) f 5' • (127) 
a=iea;jew{i) ^ Vi \d{Xi,Xj)J 

Comparing (11271) with (11031) . we obtain 

\\Cyf<M2{Dy,y), (128) 
where (7 = (7i and the constant 

2 , ,2 

M2 = — max max , 

with k{x) > Ki > 0. 

For the difference convective transport operator in the divergent form (|1 141) . 
we use the difference analogue of the Friedrichs inequality in the form (11041 ). By 
the representation (11211) . we obtain immediately: 

||C2i/f = \\Ciy + divhvyf < 2\\Ciyf + 2\\divhVy\\2. 



49 



From (11031 ) and (11041 ). for the last term in the right-hand side, we have 

||div;,t;?/||2 < ||div/,'i;||^^ ^ ^ kjd{xi,Xj) ( f l\ ] . 

Taking into account (11031) (with C = Ci), we arrive at the estimate (11281) for the 
operator C = C2, where 

2 / 

M2 = — I 2 max max + Moll div^-yll^ 



For the difference operator of convective transport in the skew- symmetric form 
(|117|) . in a similar way, we establish the inequality (11281) with the constant 



1 



M2 = — ( 3 max max \bij\'^ + Mq\\ divhv\\l^ 



Theorem 7. For the grid operators of convective transport defined in accordance 
with A114\l . dii5l) and ( \117\i . in the space of grid functions H = ^2(0;), the prop- 
erties ( lii(5D and A118\) are valid along with the estimate of the operator energy 
boundedness ( \124\i . and the estimate of subordination U28\i to the difference op- 
erator of diffusive transport U02\i hold. 

The above estimates (11241) and (11281) for difference operators of convective 
transport are fully consistent with the continuous case. They serve us as the basis 
for the study of difference convection-diffusion problems. 

3. 7. Monotone approximations on triangular grids 



Monotone approximations on triangular grids are constructed [|40|l similarly 
to other grids. We separate the positive and negative parts of the normal velocity 
component putting 

where 

Kj = lib,, + \k,\), 

Kj = ^iK -\bij\). 

To approximate the right-hand side of (II 131) . we will use the value of the grid 
function either in the central or in the peripheral node depending on the sign of 
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the velocity. This leads us to the difference operator of convective transport in the 
form 

(Cay). = - Yl hib^yj + KVi), (129) 

If we apply the difference divergence operator, then, for the difference analog 
of (|43l) . we have the expression 

* jew(i) 

Thus, we have developed the approximations (11291) and (|130l) for the convec- 
tive transport operators in the divergent (|44l) and nondivergent (l43l) forms using 
the upwind differences. 

For the boundary value problems (l88l) . (f89l) and (l89l) . (|94|) . we put into the 
correspondence the difference problems 

Cy + Dy = ip{x), x e oo, (131) 

for the grid functions y{x) = 0, x E dco. For the right-hand side of (|131l) . 
suppose, e.g., 

(p{x) = y j /(^) dx, X euj. 

To employ the fulfillment of the maximum principle at the discrete level, we 
rewrite the difference problem (11311 ) in the form 

o^iVi - = 0i, Xi E uj, (132) 

yi = 0, X E du. (133) 

Assume that ZJ is a connected grid. 

For the difference problem (I132 m T331) . the maximum principle is valid, i.e., 
the difference scheme is monotone pTI] under the following restrictions: 

ai>0, A, >0, jEWii), (134) 
5i = ai- ^ (3ij > 0, Xi E u. (135) 



51 



The difference scheme (11311 ) for the problem (|88l ). (|89l) based on the approxi- 
mations (1 1021) and (1 1301) may be written in the form (I132D . (I133D with the coeffi- 
cients ^ ^ ^ 

5,, = 0, ajj G oj. 

The monotonicity conditions (11341) . (11351) are unconditionally valid. 

Now consider the scheme (|131l) with the difference operator of convective 
transport C = Ci defined according to (11151) . The scheme (I102D . (I115D . (11311) 
may be represented in the canonical form (11321) . (11331) with 



(5j = 0, ajj G w. 
Define a local grid Peclet number as follows: 

KyXij) 

The monotonicity condition (|134l) leads to the restrictions 

Peij < 2, j G W(i), G w. (136) 

Such restrictions are typical if we apply the standard central-difference approxi- 
mations on regular grids. 

For the convection-diffusion equation with the divergent convective terms (l89l) . 
(|94|) . the use of the upwind approximations (11021) and (11291) yields 

ft, = -^'yfS + ^*(.«);^(^. J € 
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5i = divhV, Xi e u. 

Thus, the standard monotonicity conditions (11341) . (11351) are valid only if divhV > 
0. 

A similar situation occurs in the consideration of difference schemes on rect- 
angular grids. In this case, the unconditional fulfillment of the maximum prin- 
ciple for schemes with the upwind differences designed for the difference equa- 
tion (11311) may be associated with diagonal dominance by columns rather than by 
rows (as the conditions (11341) . (11351) ). The second possibility, which seems more 
promising for schemes on unstructured grids, involves the establishment of the 
maximum principle in the standard formulation for the conjugate problem. 

Consider the operator that is adjoint to C2 and is defined according to (11291) . 
Taking into account that kj = Iji, hf- = —bj^, we obtain 

Therefore ^ 

^*^'" ^ V- ^ Xi^uj. (137) 

As for the adjoint problem 

C^v + Dv = ip{x), xeu, (138) 

the unconditional fulfillment of the maximum principle written in the standard for- 
mulation is established in a usual fashion. Recall that we speak of the formulation 
for the maximum principle in the following form — if the conditions (11341) . (11351) 
are true, then the solution of the problem [1321), (11331) is non-negative (nonpositive) 
for the non-negative (nonpositive) right-hand side (11321) . 

Now we show that from the fulfillment of the maximum principle for the ad- 
joint problem, it follows that the maximum principle is satisfied for the original 
problem. For each Xi E to, we define the grid function 

f 1 

I J *^ *^ 2 5 

Shix -Xi) = < Vi 

[0, X^ Xi. 

Suppose that the grid function G{x, Xi) for the given Xi E uh the solution of the 
problem 

ClG + DG = 5h{x-Xi), XEU. (139) 
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Due to the fact that the maximum principle holds for the adjoint difference prob- 
lem (11291 ). we have G{x, Xi) > 0. 

Multiply equation (11391) scalarly by the solution of the original boundary value 
problem 

C2y + Dy = (p{x), X e 00. 
By virtu of (11391 ). the solution is represented as 

y{xi) = {G{x,Xi),if{x)). 

Thus, we get y{xi) > 0, x E to. Therefore, the maximum principle holds also for 
the original problem (fTOll) . (fT29l) . (fHTl) . 

The study of the difference scheme (11021) . (II 141) . (11311) is conducted in as sim- 
ilar way. The monotonicity of the above difference scheme (|102l) . (|1 141) . (|131l) is 
established under the restrictions (|136l) . Our investigation results in the following 
statement. 

Theorem 8. The upwind difference schemes di02P . U30\) , diiiP and di02|) . di29D . 
rtiJil) for the convection-diffusion equations U02^ , ( 17301) . ( liiiP and f |i02P . I\129\) , 
( 17371) are unconditionally monotone, whereas the schemes ( 17 021) . ( 17751) . ( I737P and 
( 17 021) . ( 17 741) . ( I737D satisfy the maximum principle under the restrictions ( 173(51) . 

The above approximations for elliptic operators of convection-diffusion are 
used for discretization in space on irregular grids for numerical solving time- 
dependent problems. 

4. Discretization in time 

Discretization in space results in the Cauchy problem for systems of ODEs 
treated as an operator-differential equation in the appropriate spaces. Two- or 
three-level difference schemes are used for numerical solving these equations. 
This part of the work discusses issues of constructing unconditionally stable schemes 
for the approximate solution of unsteady convection-diffusion problems. The in- 
vestigation is based on the general theory of stability (well-posedness) for operator- 
difference schemes. 

4.1. Two-level operator-difference schemes 

We start with the key concepts of the stability theory for operator-difference 
schemes considered in finite-dimensional Hilbert spaces. Next, for two-level dif- 
ference schemes, we formulate criteria of stability with respect to the initial data. 
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And finally, typical estimates for stability with respect to the initial data and the 
right-hand side are presented. 

For simplicity, we define a uniform grid in time as follows: 

Ur = ^r^ {T} = {t" = nT, n = 0, 1, No, tNq = T}. 

Denote hy A, B : H H linear operators in H depending, in general, on r, t". 
Consider the Cauchy problem for an operator-difference equation 

„ n+l _ „ n 

B{r)- — + = (^", e e oor, (140) 

r 

2/° = (141) 

where = yif^) E H is a desired function and (p'^,u^ E H are given. We use 
the index-free notation of the theory of difference schemes: 

y = y\ y = y''^\ y = y^-\ 

t y-y y-y 
y = , yt = ■ 

r T 
Then equation (|140|) may be written as 

Byt + Ay = ip, tE u^. (142) 

We define a two-level difference scheme as a set of the Cauchy problems 
(11401) . (11411) that depend on the parameter r. The formulation (|140l) . (11411) (as 
well as (11411) . (11421) ) is called the canonical form of two-level schemes. 

For solvability of the Cauchy problem at a new time level, it is assumed that 
B^^ does exist. Then equation (11421) may be written as 



y 



Sy + Tip, S = E -tB-^A, ip = B-'^ip, (143) 



where, as usual, E is the identity operator. The operator S is called the transition 
operator of the two-level scheme (the transition from a current time level to the 
next one). 

A two-level scheme is called stable if there exist positive constants mi and 
1712, independent of r, and cp, such that for any E H, (p E H, t E Ut, for 
the solution of (11401) . (11411) . the following estimate is valid: 

< millwl +m2 max 11^(0)11,, r E Ur, (144) 
o<e<t" 
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where || ■ || and || • ||* are some norms. The inequality (11441) reflects the continuous 
dependence of the solution of (11401) . (|141l) on the input data 
The difference scheme 

„ n+l _ „ n 

B{r)- — + = 0, re w,, (145) 

r 

l/° = M° (146) 

is called stable with respect to the initial data if for the solution of (11451) . (|146l) . 
the following estimate holds: 

<mi||M°||, eeur. (147) 

The two-level difference scheme 

„ n+l _ „ n 

B{e)^ ^ + = ip"", r e COr, (148) 

r 

y° = (149) 

is called stable with respect to the right-hand side if the solution satisfies the in- 
equality 

< ms max T G w,. (150) 

0<6»<t" 

The difference scheme (|145l) . (11461) is said to be p-stable (uniformly stable) 
with respect to the initial data in Ho if there exist constants p > and mi, inde- 
pendent of r and n, such that for any n and all ?/" G H, the solution ?/"+^ of the 
difference equation (11451 ) satisfies the estimate 

\\y"'^'\\D < pWvIId, eecor, (I5l) 

and < nil. 

In the theory of difference schemes, one of the following quantities is selected 
as p: 

P = l, 

p = 1 + CT, C > 0, 

p = exp (cr), 

where a constant c is independent of r, 

In view of (11451 ). rewrite equation (11431 ) in the form 

= Sy-. (152) 
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The requirement of p- stability is equivalent to the fulfillment of the bilateral oper- 
ator inequality 

-pD<DS < pD, (153) 

if DS is self-adjoint (DS = S*D). For an arbitrary operator of transition in (|152h . 
the condition of p-stability is given by 

S*DS < p^D. (154) 

Let us formulate the discrete analog of Gronwall's lemma. 
Lemma 5. From the estimate for the difference solution at the n + 1-st time level 

<p||2/1| + r||^"|U (155) 

it follows that the a priori estimate 

n 

<p"+i||y°||+^rp"-'^||/||, (156) 

k=0 

holds. 

Thus, from the levelwise estimate, we obtain an a priori estimate for the dif- 
ference solution at any time level. 

Let us formulate the basic criteria for stability of two-level schemes with re- 
spect to the initial data 11271 12911. The most important is the following theorem, 
proved by Samarskii, on the exact (coinciding necessary and sufficient) condition 
for stability in Ha- 

Theorem 9. Assume that in equation 17451) . the operator A is a positive and self- 
adjoint operator independing ofn. The condition 

B>-A, tecor (157) 
2 

is necessary and sufficient for stability in Ha, i-e., for the fulfillment of the estimate 

||z/"+'IU < lk°IU, teojr. (158) 

Proof. Multiplying equation (11451) scalarly by yt, we get 

{Byt,yt) + {Ay,yt) = Q. (159) 
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Using the representation 



rewrite (11591 ) as 

( (5 - yt, l/i) + ^{A{y + y),y-y)=0. (160) 

For the self-adjoint operator A, we have (Ay, y) = (y, Ay) and 

(My + y),y-y) = i^y, y) - {Ay, y). 

Substituting these relations into (11601) and using the condition (11571 ). we obtain 
the inequality 

< ||y"1U, (161) 

which ensures the desired estimate (11581) . 

To prove the necessity of the inequality (|158l) . assume that the scheme is stable 
in Ha, i.e., the inequality (11581 ) holds. We prove that this implies the operator 
inequality (11571) . Consider (11601) at the initial time level n = 0: 



2r (^(^B -^A^w,w'^ + {Ayuyi) = iAyo,yo), w 



In view of (11581) . this identity holds only if 

B - ^A^ w,u?j > 0. 

Let yo = uq ^ H he an arbitrary element, then the element w = —B~^Auo E H 
is arbitrary, too. Indeed, for any element w E H,we obtain uq = —A~^Bw E H 
since A^^ exists. Thus, the inequality holds for any w E H, i.e., we have the 
operator inequality (11571 ). □ 

The condition (11571) is necessary and sufficient for stability not only in Ha, 
but also in other norms. We now formulate (without proof) the stability result for 
Hb. 

Theorem 10. Assume that in ( \145\i . M46\i . the operators A and B are constant 
and 

B = B* >0, A = A*>0. (162) 

Then the condition di57D is necessary and sufficient for stability of the scheme 
( 17451) . I\146\) with respect to the initial data in Hb with p = I. 
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The consideration of general time-dependent problems is based on using p- 
stability. 

Theorem 11. Let A and B be constant operators and 

A = A\ B = B* >0. 

Then the condition 

^-P-B<A< i±^5 (163) 



T T 

is necessary and sufficient for the p-stability of the scheme ^7451) . U46^ in Hb, 
i.e., for the fulfilment of 

||y"+'||B<p||Z/"||B. 

Proof. Writing (|145l) in the form of (11521) . we get from (|153l) the following con- 
dition for stability in Hb'- 

-pB < B -tA< pB. 

This bilateral operator inequality can be formulated in a more traditional repre- 
sentation using inequalities in the form of (|163l) for the scheme operators. □ 

We emphasize that in this theorem there is no assumption that the operator 
A is positive (or at least non-negative). Under the additional assumption on the 
positiveness of A, we get that the condition (11631 ) is necessary and sufficient for 
the p-stability of the scheme (11451 ). (11461 ) in Ha- 

If P > 1, then stability, as in Theorem|9l is established for two-level difference 
schemes with the non-self-adjoint operator B. 

Theorem 12. Let A be a self-ajoint, positive, and constant operator. Then under 
the condition 

B > -^A, (164) 
1 + P 

the scheme ({145\) , ({146\) is p-stable in Ha- 

Proof. Adding and subtracting from the basic energy identity (see (I160|) ) 

2r ( (5 - ^ a) yt, yt) + {Ay, y) - {Ay, y) = Q (165) 



the expression 



2r^-^{Ayt,yt), 
1 + p 
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we get 



In view of (11641) and the self-adjointness of A, we obtain immediately 

{Ay, y) - p{Ay, y) + {p - l){Ay, y) < 0. 

The inequality 



with notation 



\{.Ay,y)\ < llylUlblU 
II^IU 

7] 



yields the inequality 



\y\\A 

- (p - + p < 0. 
It holds for all 1 < < p, and so we go to the desired estimate 

WvWa < WvWa, 

which ensures stability in Ha- □ 

Now we consider a priori estimates that express stability with respect to the 
right-hand side. Such estimates are employed to study convergence of difference 
schemes for time-dependent problems. 

First, we show that stability with respect to the initial data in Ho, D = 
D* > results in stability with respect to the right-hand side in the norm Hv^H* = 
\\B-'v\\r. 

Theorem 13. Assume that U40^ . ( \141\) is p-stable in Hr with respect to the initial 
data, i.e., the estimate (1757]) holds with tp^ = 0. Then the scheme A140\) . AMU is 
stable with respect to the right-hand side and the following a priori estimate is 
true: 

n 

W^^\r < P^^'Wu'Wr + 5^rp"-^||5-V1li?- (166) 

A;=0 
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Proof. Since B ^ exists, we have that equation (11401 ) may be written as 

= 5?/" + r^", S = E-tB^^A, (^" = 5-V". (167) 
From (11671 ). we get 

||l/"+1/?< ||5i/"b + r||fi-V"||ii. (168) 

The requirement of p-stability with respect to the initial data is equivalent to the 
boundedness of the norm of the transition operator S: 

||5z/"1lii<p||i/1lii, teur. 

Because of this, from (11681) . we obtain 



|l/'^+1i?^<p||l/1|/? + rl|fi-Vl 



R- 



Using the discrete analog of Gronwall's lemma, we obtain the desired estimate 
(11661) . which expresses the stability of the scheme with respect to the initial data 
and the right-hand side. □ 

In particular, li D = A ox D = B (under the condition A = A* > or 
B = B* > 0), then, from (11661) . we obtain elementary estimates for stability in 
the energy space Ha or Hb- 

Some new estimates for the two-level difference scheme (|140l) . (|141l) can be 
obtained by coarsening the stability criterion (I167I ). 

Theorem 14. Let A be a self-ajoint, positive, and constant operator and assume 
that B satisfies the condition 

B > ^^^^ (169) 

with a constant e > independing ofr. Then the scheme U40\i . AMU satisfies the 
a priori estimate 

1 -J- " 

|i,n+l||2 ^ |U,0||2 I ^ \^ ^11, „A:||2 nifW 

\\y \\a< \\u \\a + z^nv \\b-^- (170) 

k=0 
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Proof. Multiplying equation (11401) scalarly by 2Tyt, we obtain, similarly to (11651) . 
the energy identity 

2r{{B - ^-A)yt,yt) + {Ay,y) = iAy,y) + 2r{^,yt). (171) 

The right-hand side of the above expression can be estimated as 

2T{ip,yt) < 2T\\(p\\B-i\\yt\\B < 

< 2TCiiiytiii + 

with a positive constant si. Substituting this estimate into (11711) . we get 

2r yuy^+{Ay,y) < {Ay, y) + ^Ml.,. 

If the condition (1 1691) holds, then it is possible to select ei such that 

^ 1 + e, 



l-Bl 

and so 

(1 - e,)B -Ia = {1- e,){B - ii^rA) > 0, 

{Ay,y) < (A2/,2/) + i^r||y;|||_i. 

The last inequality implies the estimate (11701) . □ 

Theorem 15. Let Abe a self-ajoint, positive, and constant operator, and assume 
that B satisfies the condition 

B>G+^A, G = G* >0. (172) 

Then the solution of U40\i . AMU satisfies the a priori estimate 

n 

•I ,n+l||2 ^ |U,0||2 I \^^||,„A:||2 /'1'70\ 

\y lU < If lU + 2 Z^^llv' IIg-1- (1'3) 



k=0 



62 



Proof. In the identity (11711) . we employ the estimate 



Substituting this estimate into (11711) and taking into account (11721) . we get 

{Ay,y) < {Ay,y) + K\\^\\l-^ 

that, by a discrete analog of Gronwall's lemma, gives (|173l) . □ 

The convergence study of difference schemes is conducted in various classes 
of smoothness of the solution of the original differential problem, and therefore we 
must have a wide range of estimates. In particular, the right-hand side should be 
estimated in different and simply calculated norms. Only typical a priori estimates 
for solutions of operator-difference schemes are considered here. 

We now apply the above results to elementary schemes with weights for an 
operator-differential equation of first order. The Cauchy problem 

— + Au = f{t), t>0, (174) 

m(0) = u\ (175) 
with A >0 is associated with the two-level scheme with weights 

^ ^ + A(ay"+l + (l-a)l/") = y.^ reur, (176) 

T 

y° = (177) 

The scheme (11761 ), (11771 ) may be written in the canonical form (11421) with the 
operators 

B = E + arA, A > 0. (178) 

Theorem 16. The scheme with weights ( \176\i . ( \177\i is stable in H with respect 
to the initial data if and only if the following operator inequality holds: 

A* + (a -]AtA*A>Q. (179) 
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Proof. By A > 0, there exists A~'^. Multipying (11761 ) by we go from (|142l) . 
(11781 ) to the scheme 

— + iy" = r G UJr 

T 

where 

5 = + arE, A = E. 

The necessary and sufficient condition for stability of this scheme with respect to 
the initial data iwH = (Theorem IH) is formulated as the inequality 

A^^ + -^tE> 0. 

Multiplying it from the left by A* and from the right by A, we obtain (I179I ). This 
completes the proof of the theorem. □ 

If a > 0.5, then the operator-difference scheme (11761) . (11771) is unconditionally 
stable (stable for any r > 0). 

4.2. Difference schemes for convection-diffusion poblems 

Discretization in space of the Cauchy problem (|37|) . (HTI) yields the problem 
(see, e.g., dlB, dSll)): 

^ + Ay = ^{t), A = C + D, 0<t<T, (180) 
at 

y{0) = u^. (181) 

With the above approximations, the grid operators of convective and diffusive 
transport inherit the basic properties of differential operators in the appropriate 
spaces of grid functions. Among these properties, we recall the following as 
the major ones. The constant (time-independent) grid diffusion operator is self- 
adjoint and positive definite: 

^D = D^ D = D*>^K,E, K,>0, Mo>0. (182) 
at at Mq 

For the grid operator of convective transport in various forms (C = C{t) = 
Ca, a = 0, 1, 2), we have 

Co = -C;, (183) 
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\{C^y,y)\<M4y\\\ « = 1,2, (184) 

\\C^yf<M2{Dy,y), a = 0,1, 2, (185) 

with the corresponding positive constants Mi and M2. 

To solve numerically the problem (11801) . (|181l) . we consider the two-level 
scheme with weights 



n+l _ 



y ^ -y . .n+i 



(186) 

+ Z^(a2i/"+i + (1 - a2)y") = </^^ r e a;., 

yo = Uq. (187) 

Here, e.g., we have 

C = C(0.5(r+^ + r)), yp" = (^(0.5(r+^ + r)). 

Among the most important variants of the difference scheme with weights 
(11861) . (11871) . we highlight the scheme with equal weights (cti = (T2) and the 
scheme, where convective transport is taken from a previous time level (cti = 0). 

We start with the convective transport operator in the skew- symmetric form, 
i.e., C = —C* = Cq. Problems with the convective transport operator in the 
nondivergent (C = Ci) and divergent (C = C2) forms will be considered later. 
Assume that in the difference scheme (|186|) . we have 

(Jl = CT2 = (T. (188) 

In view of (11881) . instead of (11861) . we consider the difference scheme 

^ ^ + (Co + + (1 - a)y^) = yp", e e ujr. (189) 

r 

The scheme (|187l) . (11891) under investigation may be written in the canonical 
form for the two-level difference scheme (11401) . (11411) with the operators 

B = E + arA, A = Co + D>0. (190) 

The main peculiarity of difference schemes for the convection-diffusion equation 
is connected with non-self-adjointness of the operators B and A. Therefore, it is 
impossible to use the above results on stability of operator-difference schemes, 
which were formulated for constant self-adjoint operators. 
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The second important feature is associated with the fact that operators of the 
difference scheme are variable in time. We consider the problems with the time- 
dependent difference operator of convective transport. To obtain a priori estimates 
for such problems, it is often necessary to require additionally Lipschitz continuity 
of the difference operators with respect to time. 

Conditions for stability of the scheme (11401 ). (11411) . (11901 ) have been presented 
above in the form of Theorem \T6\ Let us supplement this result with the corre- 
sponding stability estimate of the difference solution with respect to the right-hand 
side and the initial data. 

Theorem 17. The difference scheme di40D . A141\) . A190\) is unconditionally stable 
for a > 0.5, and the difference solution satisfies the a priori estimate 

1 " 

iiy"+^ir<ii^°f (191) 

fc=0 

Proof. Rewrite the scheme (I140L (1 1901) (see (11891) ) as follows: 



yn+l _ yu 
T 



At;"+i = (^", n = 0,l,...,iVo-l, (192) 



where 



= + (1 - ay") ={a-l]yt + ^(2/"+^ + 



yTl+l _ yU 

yt = • 

r 

Multiplying equation (11921) scalarly by f "+^, we obtain 
r{yt,yt) + {Av-^\v-+') 

ZT 



(193) 



In the condition (11901) . we have (Ay, y) = (Dy, y). For the right-hand side, we 
use the estimate 
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With this in mind, from (11931 ). under the conditions of the theorem we get the 
estimate 

Thus, we come to the desired estimate (11911) . □ 

The a priori estimate (11911) obtained above for the difference solution is a grid 
analog of the a priori estimate (l60l) for the solution of the differential problem 
(|4TI) . (l59l) . because the convective transport operator in the skew- symmetric form 
under consideration corresponds to the constant A^i = in (l60l) . 

Now we consider the case, where the skew-symmetry of the difference op- 
erator of convective transport is not valid. We will study the problem with the 
convective transport written in the nondivergent form, i.e., C = Ci. The case the 
convective transport in the divergence form (C = C2) is investigated in a similar 
way. 

Let us examine the scheme (11401) . (11411) . where 

B = E + arA, A = Ci + D. (194) 

It is important to distinguish two classes of problems. The simplest case is asso- 
ciated with the assumption that the operator A is non-negative. Such a situation 
takes place, e.g., if MiMq — ki < — convective transport has only an insignifi- 
cant effect. Indeed, in view of (11821) . (11841) . in the case (11941 ). we have 

{Ay,y) = {Ciy,y) + {Dy,y) > -Mi\\yf + T^«:i||l/f 

IVIq 

= ^K-MiMo)||?/|r. 
Because of this, for the operator A, we have the following lower bound: 

A>^{ki-M,Mo)E. (195) 

Mq 

Another case deals with slightly compressible flows, where A > under the 
condition M2M0 — ki < 0. In this situation, in view of (11821 ) and (11851 ). we obtain 

{Ay,y) = {Ciy,y) + {Dy,y) > ||y|| + {Dy,y) 
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Under these restrictions on parameters of the problem, we can apply the results on 
the unconditional stability (Theorem [161) for the difference scheme (11401 ). (|141l) . 
(fT94l) in Hfora> 0.5. 

In the general case, we cannot rely on the non-negativity of the operator A. 
This leads to the fact that the conventional schemes with weights are not uncon- 
ditionally stable under the standard restrictions a > 0.5. Let us consider the 
difference scheme (11401) , (11411) . (11941 ) as an example. 

The solvability of the scheme (11401) . (11411) . (11941) (B > 0), in view of the 
fact that the operator A is not non-negative, takes place under the constraint of an 
appropriately small time step — we speak of conditional solvability. Taking into 
account (11941 ). (11951) with MiMq — ki > 0, we get the following restriction on a 
time step: 

r<n = -. (196) 

(t(MiMo-k) 

In this case (see Theorem [21 the estimate ([60l) for the solution of the differential 
problem), it is necessary to be oriented to obtaining an appropriate estimate that 
expresses conditions for f)-stability. 

We have already formulated the necessary and sufficient condition for ^-stability 
in the case with the constant self-adjoint operators B and A. Therefore, our study 
will be based on the schemes with weights of type (11401) . (11411) . (11941 ) considered 
above. 

Let us define new grid functions f 

y" = f,V, n = 0,l,...,iVo, g>0. (197) 

A condition for f)-stability for is evidently equivalent to stability (g = 1) for 
f Substitution of (11971) into (11451) yields the difference scheme 

B + Av"" = 0, tnE COr, (198) 

T 

where 

B = qE + agTA, A = + {1 + a{g - 1))A. (199) 

r 

It is possible to use the following representation for the operators of the difference 
scheme (11981) : 

B = G + arA, (200) 
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which treat the scheme (11981) as a scheme with weights. In view of the represen- 
tation (11991) . we obtain in (I200D : 



G = y -E, d = ^ -. (201) 

l + (x{g-l) ' l + 

Similarly to Theorem [161 we prove the stability of the scheme (11981 ). (|200l) in 
Hg, i.e., in H with a > 0.5 under the constraint A > 0. Taking into account 
(12011 ). we get the desired condition on a weight of the difference scheme (11981 ). 
(fT99l) : 

a > (202) 

1 + g 

The non-negativity of the operator A is connected with an appropriate choice 
of g. In view of the stability estimate for the differential problem (see Theorem|2l 
the estimate (l60l)). it is natural to set 

g = l + MiT. (203) 

Taking into account the estimate (|195l) . the condition A> (see (11991) ) is fuUfiled 
for 

Ml Mo - (1 + (xtMi){MiMo - /ti) > 0. 
This inequality yields the following restriction on a permissible time step: 

^ - = aM,iM?Mo-n,y ^^^^^ 

A comparison with the estimate (11961) shows that the time step restriction (12041) is 
slightly stronger {t2 < ti, we recall, MiMq > /ti). Summarizing, we obtain the 
following statement. 

Theorem 18. The scheme with weights di40l) . {\141^ . U94^ under the constraint 
MiMq > Ki is g-stable in H, where g is defined according to ( I203D . if the weight 
a satisfies the restriction Ii202\) and a time step meets the condition A204\) . 

This statement complements Theorem [l7l which ensures the stability of the 
scheme (11401) . (|14H) . (|194l) under the constraint MiMq < ki in H with a > 0.5. 
Possible non-negativity of the operator A = Ci + D leads to the situation, where 
we must use f)-stability. In addition, we impose (see (|204l) ) restrictions on a time 
step. 
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In solving convection-diffusion problems, it is reasonable to focus on differ- 
ence schemes, where a part of the operator A (it is, of course, the convective trans- 
port operator) is taken from the previous time level [|32n . Such explicit-implicit 
schemes from the above class of two-level schemes with weights are considered 
in n3% . Suppose now that in the difference scheme (I186|) . we have 

(Xi = 0, (T2 = a. (205) 

The homogeneous ((p„ = 0) scheme (11861 ). (12051) is reduced to the canonical form 
(11401) if we define 

B = E + arD, A = C + D. (206) 

For any r > 0, we have B > 0, and therefore the discrete equation (11861) . (12051) is 
solvable at every time level. Let us formulate a sufficient condition for f)-stability 
of the difference scheme for the convection-diffusion equation in Ho. 

Theorem 19. The solution of the explicit-implicit scheme U86^ , ( 12051) with a > 
0.5 satisfies the estimate 

||l/"+iD < qWWd (207) 

where 

g = l + -fr, (208) 
and M2 is the constant from the inequality A185\l . 

Proof Multiply (11861) scalarly by 2ryt = 2(y"+^ - y") and, in view of (|2061) . 
obtain the energy identity 

r{{2B - TD)yt, y^) + (Z}y"+\ - (Dy", y") + 2r(Cy", y^) = 0. (209) 

Taking into account the representation (|206l) and the constraint a > 0.5, from 
(12091 ). it follows the inequality 

2T{y„y,) + {Dy-+\y-+') - {Dy\y-) < 2r\{Cy- ,y,)\. (210) 

In view of (11851 ). the right-hand side is evaluated as follows: 

\{Cy-.yt)\ < btf + IwCy'-r < btf + ^ 
Substitution into (12101 ) yields 
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Therefore, in view of inequality 

M2 / M2 ^ ^ 

we obtain the desired stability estimate (12071) . (I208I ). □ 

The f)-stability estimate (|207l) . (12081) . derived here, is fully consistent with the 
corresponding estimate for the differential problem (see, e.g., the estimate (|62l) 
and the proof of Theorem [2l). An important point is that, in contrast to Theo- 
rem [T8l we obtained stability with the standard restrictions on a weight cr in a 
stronger norm. Moreover, the implementation of the explicit-implicit scheme is 
much simpler from the computational point of view - we must invert a self-adjoint 
elliptic grid operator. 

Considering two-level difference schemes, we have highlighted two main classes 
of difference schemes for unsteady convection-diffusion problems. The first class 
is based on the use of the simplest schemes with equal weights for the convective 
and diffusive transport. The second and the most promising class of difference 
schemes (explicit-implicit schemes) is associated with the explicit treatment of 
the convective transport. Here we do not analyze the complete set of three-level 
difference schemes. We focus on the study of explicit-implicit schemes. Using 
three-level difference schemes, we can obtain the second-order approximation in 
time. 

To solve numerically the problem (|180l) . (|181l) . we employ the three-level 
explicit-implicit scheme with weights 

„ n+l _ n,n~l 

^ ^/ + Diay-^' + (1 - 2a)y- + ay^-') ^^^^^ 
+ Cy- = ^\ n=l,2,...,iVo-l 

with 

^ = ^0, y^=u\ (212) 

In (|211l) . we put, e.g., C = C{t"), v?" = v'(t"). To specify the second initial 
condition {u^ in (12121) ) with the second order, in the simplest case, we involve a 
two-level scheme, so that 

^ + (C + Z^)^ = /. 

The difference scheme (121 II) . (12121) approximates (11801 ). (11811 ) with the second 
order in time. 
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The explicit-implicit scheme (121 II) is written in the canonical form 

B{r f + - 2y- + y--') ^^^^^ 

+ = <^", n = l,2,...,No-l 

with 

B = E, R = aD, A = C + D. (214) 

To evaluate the difference solution, we introduce the norm associated only with 
the diffusive transport, i.e., 

( W (215) 

Stability is established taking into account the subordination of the convective 
transport operator to the diffusive transport operator — we say about the estimate 
(fT85]) . 

Theorem 20. If o > 0.25, then the difference scheme ( 127 7 D . ( \212\) is g- stable with 

g=l + M2-^T, (216) 
4cr — 1 

and the solution satisfies the a priori estimate 

^"+1 < ^f" + r||^"||2. (217) 
Proof. For the scheme (12131 ). (12141) . we have 

2r 

where 

For the first two terms in the right-hand side, we obtain 

|(Ci/",u7"+^ + w")| < -^||u;"+i+w"f + r||Cy"f, 

4r 
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Thus, in view of (12151) . we arrive at the inequality 

gn+i < + rMally"!!^ + r||(^"f . (218) 
Next, we use the estimate 

An 

IIz/"IId < (219) 

- 4cT- 1 

According to (12121) . we get 



= {Dy^, y-) - r(Dy", y,) + ar^Dy,, y,) 

|2 

Id 



>{^-myl\l + {^-^y\\yi\\l. 

For a > 0.25, we select the parameter (3 = 1/ (4cr) and obtain the estimate (|219l) . 
Substitution of (12191) into (12181) yields the levelwise estimate (12161) . (12181) . □ 

4.3. Unconditionally stable schemes 

For convection-diffusion problems with convective transport in the divergent 
and nondivergent forms, we have constructed (Theorem [T8l) conditionally stable 
schemes with weights. Restrictions on a time step (see ( 12041 )) are governed by fea- 
tures of the problem and do not related, in general, with parameters of discretiza- 
tion in space. Conditionally stable schemes with weights are developed only for 
problems with the convective transport in the skew-symmetric form (Theorem 

Different nature of convective and diffusive transport as well as reaction pro- 
cesses appear, in particular, in significantly distinct representative rates of these 
phenomena. Such heterogeneity can be taken into account when choosing dis- 
cretization in time. The most pronounced occurrence of the heterogeneity of 
discretization in time is expressed in explicit-implicit schemes. In this case, for 
numerical solving the unsteady problem, a part of the problem operator terms is 
approximated by explicit relationships, whereas the other part is treated implicitly. 

Explicit-implicit schemes are widely used for the numerical solution of convection- 
diffusion problems. Various variants of inhomogeneous discretization in time are 
given in djl. One or another explicit approximations are applied to the convec- 
tive transport operator, whereas the diffusive transport operator is approximated 
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implicitly. Thus, the most severe restrictions on a time step due to diffusion are 
removed. In view of the subordination of the convective transport operator to 
the diffusive transport operator, we have already proved unconditional stability 
of the above-considered explicit-implicit schemes for time-dependent convection- 
diffusion problems. 

Similar techniques are used in the analysis of diffusion-reaction problems. In 
this case (see, e.g., [|24l1 ). the diffusive transport is treated implicitly, whereas 
for reactions (source terms), explicit approximations are used. Such explicit ap- 
proximations demonstrate obvious advantages for problems with nonlinear terms 
describing reaction processes. 

In convection-diffusion-reaction problems, the problem operator may be sign- 
indefinite. This means that the system may be nondissipative, i.e., the solution 
norm for the homogeneous problem does not decrease during the time evolution. 
Thus, the exponential growth of the solution may be observed, and such behav- 
ior of the solution must be reflected at the discrete level. Unconditionally stable 
schemes for such problems are constructed in the work [41] . They are based on 
the splitting of the problem operator into two terms, where one of the terms has 
explicit approximations in time, whereas the other is approximated implicitly. Im- 
plicit approximations are applied to the part of the problem operator that causes 
the dissipative properties of the problem. In the case of the skew-symmetric oper- 
ator of convective transport, such a splitting is used for the operator of reaction. 

The standard schemes, which are used in computational practice, should be 
corrected even for solving dissipative problems. For example, both the standard 
fully implicit scheme (backward Euler) and symmetric scheme (Crank- Nicholson) 
does not produce the exact solution for the test problem (A > 0): 

^ + Am = 0, u(0) = M°. 
at 

In ifnll . there is discussed a modification of standard schemes that is based on the 
use of (exp(Ar) — 1)/A instead of the original time step r in the application to 
the fully implicit scheme. More recent results concerned with constructing and 



employing such nonstandard discretizations in time can be found, e.g., in [11811 . 
Here we mention new possibilities in designing unconditionally stable schemes 
for solving unsteady convection-diffusion problems that involve the introduction 
of new variables. 

Time-dependent convection-diffusion problems with the convective transport 
in the divergent (I37l)-(l39l) and the nondivergent (I36l)-(l38l) forms may be written 
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as the Cauchy problem for the operator equation (compare with (|4T1) ): 

+ Au = f{t), A = Co + n + V. (220) 

at 

Here we introduce the reaction operator 

TZu = r{x, t)u. 

In the case dSTJ-dlll), we have 

r{x, t) = —- div V. 



Similarly, for equation (I36l)-(l38l). we obtain 

r{x, t) = - diyv. 
For the reaction operator, we get 

n = n\ ms <n< ms. (221) 

Using roughened estimates fr the reaction operator, we can put 

m = -Ml, M = Ml. 
After discretization in space, from (12201) . we obtain the equation 

^ + Ay = ip{t), A = Co + R + D, 0<t<T, (222) 

supplemented by the initial condition (|181l) .. 
For the operator R, we have 

Ry = r{x,t)y, xEu. (223) 

In this case, we get 

R = R\ mE <R< ME. (224) 

For instance, the convective transport operator of in the nondivergent form seems 
like this: ^ 

r{x,t) = --divhV 
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for an appropriate approximation of the divergence operator. 

To construct unconditionally stable schemes for solving the problem (11811) . 
(12221) without the assumption of non-negativity of the problem operator, we apply 
explicit-implicit approximations. The bottleneck is connected with the reaction 
operator, and therefore for m < 0, we split it into two terms: 

R = R+ + R_, R+ = Rl, R_ =R*_, < R+ < ME, mE < R^ < 0. 

(225) 

By (12231 ), it is sufficient to put 

R+y = f+i^x, t)y, R^y = r_(cc, t)y, x eu, 

where 

r+ = max(0, r), r = r+ + r_. 

Using two-level explicit-implicit schemes, we may rely only on the first-order 
accuracy with respect to time. Therefore, we focus on the fully implicit approxi- 
mations of the main operator terms. We employ the difference scheme 

— + (C" + + RDy""-^^ + Rly"" = 0, n = 0, 1, A^o - 1- (226) 

T 

Theorem 21. The explicit-implicit scheme U41\) . ( \223\) -( ^226\) with m < is un- 
conditionally g- stable in H for 

g = 1 — mr, (227) 

and the difference solution satisfies the estimate 

< n = 0,l,...,iVo-l. (228) 

Proof. Multiplying equation (12261) scalarly in H by y""*"^, and taking into account 
the skew-symmetry of the operator Co, positive definiteness of the operator D, 
and relation (12251) . we obtain 

||y"+i'<(2/"+\2/'^)-(W,2/"). (229) 

In view of 

\-{Rly\yn\<My''^'\m 
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from (12291 ). we arrive at the inequality 

||^n+i||2 < (1 _2mr)||y"f . 

By virtu of 

(1 - 2mT) < (1 - mrf, 
this yields the inequality of £»- stability (12281) with g defined by (12271) . □ 

Among possible generalizations of explicit-implicit scheme (I141|) . (12261) . spe- 
cial attention should be given to schemes of the second-order accuracy with re- 
spect to time. The symmetric scheme provides an example of such a scheme: 

„ n+l _ „,n~l n,n+l _i_ n n < „,n+l 

+ {C +D + R^) ^230) 

+ i?V = 0, n = l,2,...,iVo-l, 

where now C" = C()f:"), i?" = R{t"). To start calculations with the second order 
in time, we put, e.g., 

+ liiC + D + R')y' + (C'' + D + RV) = 0, 
r 2 

Because of this and taking into account the initial condition (11411) . the difference 
equation (12301) is considered for a given and y^. 

In addition to (12301 ). special mention should be given to the scheme 

+ Rl{2f' - = 0, n = 1, 2, iVo - 1. 

Preserving the second-order approximation in time, for this scheme, the implicity 
of the main part of the problem operator is expressed more essentially. 
In equation (12221) . for the operator A, by (|182l) and (|183|) . we have 

A > mE + —KiE. 
Mo 

In our study, we use a more rough estimate 

A > mE, (231) 
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and consider the most interesting case m < 0. 

To construct unconditionally stable schemes for the differential problem (|141l) . 
(I222|) . (|224l) under the condition (12311) . we define a new function w: 

y = exp{—mt)w. (232) 

Substitution of (|232l) into (11411) . (12221) for the homogeneous right-hand side gives 
the following problem for w: 

^ + Aw = 0, A = A- mE, < t < T, (233) 
dt 

w (0) = (234) 

For this transformation, the problem operator A is non-negative (A > 0). 

To solve the problem (12331) . (12341 ). we apply the two-level scheme with weights: 



+ i"((Tu;"+i + (l-(T)u;") =0, recor, (235) 



w° = (236) 

This scheme under the standard constraints a > 0.5 is unconditionally stable 
(Theorem [TtI). 

Let us write the difference equation (I235|) for the desired grid function y"-. 
Taking into account = + r, we put 

= exp(— mt"')^"', y""^^ = exp(— mt") exp(— mr)ty""^^ 

Because of this, from (12351 ). (12361) . we obtain the following difference scheme for 

exp(mr)7/"+i - ^ _ (^exp(mr)y"+i + (1 - a) y") = 0, (237) 

T 

y° = (238) 



In contrast to the nonstandard schemes discussed in [Il7l.ll8ll. a positive effect 
is achieved not only through the use of new approximations in time, but also by 
correcting the problem operator. 
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Theorem 22. The difference scheme t\237\) , ({238\) for a > 0.5 is unconditionally 
g-stable in the H with 

Q = exp(— mr), (239) 
and the solution satisfies the a priori estimate 

\\y'''-'\\<Q\\yl\. (240) 

Proof. The above proofs were based on the transition to the problem with a non- 
negative operator and the use of the previous Theorem[T7l It is possible to conduct 
a direct proof of stability for the scheme (12371) . (12381) . Rewrite the scheme under 
consideration in the form 

exp(mr)y"+^ — , ^ n+i 



T 



+ = 0, re Ur, (241) 



where 



= aexp(mr)y"+i + (1 - a) 



exp(mr)?/""'""'^ — 



r 

Multiplying equation (|241l) scalarly by p""*"^, we obtain 

+ ^ ((exp(mr)i/"+\exp(mr)i/"+i) - = 0. 

From this equation, under the conditions a > 0.5 and A > 0, it follows that the 
stability estimate (I239L (12401) holds. □ 

It is important to note that, in contrast to the conventional scheme with weights 
(see Theorem[T8l), here stability is obtained with no restriction on a time step. The 
value of g defined by (12391 ) is fully consistent with the corresponding constant for 
the solution of the differential problem. The transition to a new time level involves 
the solution of the grid problem 

{E + aT{A - mE))y''+^ = x". (242) 

The equation (12421) is a system of linear algebraic equations with a positive definite 
and non-self-adjoint matrix; it can be solved using standard iterative methods. 
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5. Stability in Banach spaces 

The main results on stability of difference schemes for the unsteady convection- 
diffusion equation were obtained above considering the problem in Hilbert spaces 
of grid functions. Here we study difference schemes in Banach spaces, where 
stability of difference schemes is established in the uniform and integral norms. 

In our study we can employ the maximum principle for difference schemes 
as it was done in investigating monotone approximations. The second and more 
promising approach presented below is to use the concept of the logarithmic norm. 
In this section, monotone schemes of the second-order accuracy in space are con- 
structed for the time-dependent convection-diffusion. 

5.1. One-dimensional problems 

To simplify the material presented here, we start with the ID convection- 
diffusion problems. Consider the time-dependent convection-diffusion equation 
with convective terms in the nondivergent form: 

for 

0<x</, 0<t<T. 
This equation is supplemented with homogeneous Dirichlet boundary conditions: 

m(0, t) = 0, u{l, t) = 0, 0<t<T. (244) 

In addition, the initial condition is given: 

u{x, 0) = u\x), 0<x <l. (245) 

The second important example is the unsteady equation of convection-diffusion 
in the divergent form: 

Consider the set of functions u{x, t) satisfying the boundary conditions (|244l) . 
The transient problem of convection-diffusion is written in the form of the operator- 
differential equation 

+ Au = f{t), A = A{t) =C{t) +V, (247) 
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where C{t) denotes the convective transport operator, and V stands for the op- 
erator of diffusive transport. The Cauchy problem for the evolutionary equation 
(12471) is supplemented with the initial condition 

u{0) = (248) 

We recall some a priori estimates for the convection-diffusion problems (12431) - 
(12451) and (|244I) - (I246I) . which are derived from the maximum principle. The corre- 
sponding a priori estimates are derived in the spaces Loo(0, 1) and Li(0, 1), where 
the norms are, respectively, 

||f lloo = max |'y(a:)|, \\v\\i = / \v{x)\dx. 

0<x<l Jq 

The solution of the time-dependent convection-diffusion problem (I243l) - (|245l) 
(the convective transport in the nondivergent form) satisfies the a priori estimate 
in Loo (0,/): 

Mx,t)\\^<\\u\x)\\^+ [ \\f{x,e)\\^de. (249) 

We present also the estimate for the convection-diffusion equation with convective 
terms in the divergent form. The solution of the problem (I244I) - (I246I) satisfies the 
a priori estimate in Li(0, /): 

lk(x,t)||i< ||n°(x)||i+ / \\f{x,e)\\rde. (250) 

The a priori estimates (12491) , (12501) serve us as a guide in considering discrete 
problems. 

5.2. Stability of two-level schemes 

Let us obtain sufficient conditions for the stability of two-level difference 
schemes for the Cauchy problem for a system of ODEs. Further, these general 
conditions will be applied to particular cases of model convection-diffusion equa- 
tions with the convective terms in the nondivergent and divergent forms. 

Consider a system of linear ODEs of first order: 

d 

-j^ + '^aij(t)wj = (piit), i = l,2,...,m. (251) 
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Assume that w = w{t) = {wi, W2-, • • • , Wm}-, A = [aij], then we can write (12511) 
in matrix (operator) form as 

^ + A{t)w = 4>{t). (252) 
at 

We will construct difference schemes for numerical solving the Cauchy problem 
(12521) for t > and the initial condition 

w{0) = (253) 

We will investigate the stability of the difference solution of the problem (|252l) . 
(12531) in Loo and Li. For a norm of a vector and a norm of a matrix, consistent 
with it in L^, we have 



\w\ 



max \wi\, \\A\\oo = max N \aij\. (254) 

l<i<m l<i<m ' ^ 



Similarly, in Li, we obtain 

m m 

ll^lli — / ll^lli = max y |ajJ. (255) 

1=1 j=l 

The problem (12521 ), (12531) will be considered under the following constraints. 
Assume that the diagonal elements of the matrix A are non-negative, and there is 
row-wise or column-wise diagonal dominance, i.e., we have 

flii > ^ Iflijl, 2 = 1,2, ...,m (256) 

(weak diagonal dominance by rows) or 

On > ^ Ictjil, z = l,2, ...,m (257) 

(weak diagonal dominance by columns). 

The logarithmic norm of the matrix A is defined |[8i[i2] by the number 

ii\A\ = Mm ^f^^^ -. 

5^0+ 5 
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For the logarithmic norm of a matrix in L^o (consistent with (12541 )) and in Li 
(consistent with (12551 )). we have the expressions 



m 

/ioo [A] = max (aii+ 

^ ' l<i<m \ ^-^ 



m 

\a 



m 

lii[A] = max (ajj + ^ 

l<7<m \ ' 



m 

'a,; 



l<j<m 

In view of the restrictions (12561) . (|257l) . we have that the logarithmic norm of the 
matrix —A in the Cauchy problem (12521) . (|253l) satisfies the inequality 

A-A\ < (258) 

in the corresponding space (in for (|256l) and in L\ for (12571) ). 

Among the properties of the logarithmic norm (see iSl, |9|]), we highlight the 
following: 

1. /i[c74] = c/i[y4], c = const > 0; 

2. \i\cE + y4] = C + yU[A], C = COHSt; 

3. \Aw\ > max{— y4], — yu[y4]} ||w||. 

The emphasis is placed on the property 3, which allows to get easily the lower 
bound of the norm Aw. This bound can be combined with the standard upper 
bound of Aw: \\Aw\\ < \\A\\ \\w\\. 

Let us study the stability of difference schemes for the problem (12521) . (|253l) . 
We denote the approximate solution at the time level = nr (where r is a time 
step) as y", and write the two-level difference scheme with weights 



„ ra+l _ „,n 

^ ^ +A{ay^+' + {l-a)y^)=<f^, (259) 



T 



where, e.g., A = A(crt"+^ + (1 - 0")^"), with the initial data 

y° = u^. (260) 

A sufficient condition for stability of the scheme (12591) . (12601) is formulated as the 
following statement. 
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Theorem 23. Assume that in the Cauchy problem ^252P . ( I253P the matrix A sat- 
isfies the restriction ^25(5P (or ^2571) ). Then the difference scheme with weights 
A259\) . ^2(501) z^' unconditionally stable for a = 1, and it is conditionally stable for 
a < 1 in LoD (in Li) if and only if 

T < ( max ttii] . (261) 



1 — a V i<j< 

In this case, the difference solution satisfies the a priori estimate 

n 

< ||n°|| +^r||(^'=||. (262) 

fc=0 

Proof From (1259k it follows that 

{E + arA)y"+i = {E-{1- a)TA)y'' + ry?", 

and therefore 

\\{E + aTA)y"+^\\ < ||(^ - (1 - a)rA)?/"|| + r||(^"||. (263) 

For the left-hand side of (|263l) . by the above-mentioned properties of the logarith- 
mic norm and in view of (|258l) . we have 

\\{E + arA)?/"+i > -^i[~E - arA] 

= {l + a^,[-A])\\y-+'\\>U^\ 

For the first term in the right-hand side of (|263l) . we obtain 

Wi^E - (1 - a)TA)y-\\ < \\E - (1 - a)TA\\ m 

We investigate this estimate in more detail for Loo- The case Li is studied in 
a similar manner. Considering (|254l) and taking into account the condition of 
diagonal dominance (|256l) ). we have 

m 



a,: 



\E — (1 — a)TA\\ = max 1 — (1 — cr)r ( ajj + \ 

l<i<m \ ^ 

< max ( |1 — (1 — a)Taii \ + (1 — a)T N 

l<i<m \ ' 



m 



< max (|1 - (1 - cr)Taii\ + (1 - (T)rajj) < 1 

l<i<m 
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with < (T < 1 and under the restriction (12611 ) on the time step. 
The substitution into (12631 ) yields the inequality 

< ||y"1| + r||^"||, 

which immediately implies the desired estimate (12621 ) for stability with respect to 
the right-hand side and the initial data. □ 

The above estimates for stability (12621) in Loo and Li are directly associated 
with the monotonicity of the difference solution of the problem [2591) . (12601) under 
the assumption that the off-diagonal elements of the matrix A are non-positive. 
Let us prove the following statement. 

Theorem 24. Assume that in the schemes ( 12591) . ^260i . the conditions of diagonal 
dominance ( 125(51) (or ( I257D ) are fulfilled for 

aij < 0, j, i, j = 1, 2, . . . , m (264) 

and let 

u°>0, <^">0, n = 0,l,..., 

then 

y"+i>0, n = l,2,..., 

for any t > if a = 1, and z/0 < a < 1, this is true under the constraints on a 
time step A261\) . 

Proof. For the transition from the current time level to the next one, we have 

+ aTAy"+' = g^, n = 0, 1, . . . , (265) 

where 

g^ = y^-(l- a)TAy'' + Tip''. (266) 

Suppose that > (for n = this is true from the assumptions of the theorem). 
We show that from this it follows also the non-negativity of y""*"^ (y"'*'^ > 0). 

We prove that under the assumptions of the diagonal dominance (12561 ) (or 
(12571) ) and under the restrictions on a time step (12611) . for a non-negative and 
93", we get (7" > 0. . In view of (|266l) . we obtain 

m 

= {1 - (1 - a)Ta,,)y^ - [1 - cr)r ^^jV? 
> (1 - (1 - fx)ra,,)yr > 0. 
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In the conditions of the theorem, the matrix of the system of linear algebraic 
equations (12651) is an M-matrix, i.e., we have: strong diagonal dominance, positive 
diagonal elements, and non-positive off-diagonal elements of the matrix. Because 
of this, from g"" > 0, it follows that ?/"+^ > 0. □ 

Apply the derived results to studying stability and monotonicity of difference 
schemes for time-dependent problems of convection-diffusion in the nondivergent 
and divergent forms. 

5.3. Difference schemes for convection-diffusion equations 

For simplicity, we restrict ourselves to uniform grids. On the interval [0, /], we 
introduce a grid 

u! = u U du = {x \ X = Xi = ih, i = 0,1, . . . , N, Nh = I}, 

where to is the set of interior nodes: 

u = {x \ X = Xi = ih, i = 1,2, . . . , N — 1, Nh = I}. 

After discretization in space of the model convection-diffusion problems with 
homogeneous boundary conditions (I243I) - (I245I) and (I244I) - (I246I) . we arrive at the 
problem (12521) . (12531) . where m = N — 1 and the approximate solution Wi{t) = 
w{x, t), X E CO. The difference diffusion operator is specified, e.g., as follows: 

Dw = — -r^k{x + 0.5h){w{x + h,t) — w{x, t)) 

'J (267) 
+ T-^-fcfx — 0.5h)(w(x, t) — w(x — h,t)), X E u 

with 

w{x,t) = 0, xedu. (268) 

Approximation of convective transport is conducted in such a way that v{x, t) 
are defined at the half-integer grid points u. For operators of convective transport 
in the nondivergent form (equation (|243l) ). in view of (|244l) . we put 

Cw = -^v{x + 0.5h, t){w{x + h,t) — w{x, t)) 

'^^ (269) 

-I — -v(x — 0.5h, t)(w(x, t) — w(x — h,t)), x G u. 
2h 
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A similar approximation of the second order with respect to h for the convective 
transport operator in the nondivergent form (equation (12461) ) leads to 

Cw = -^Tj-v^x + 0.5h, t){w{x + h,t) + w{x,t)) 

(270) 

— t7tv{x — O.bh, t){w{x, t) + w{x — h,t)), x E u. 

Let us formulate the condition for stability and monotonicity of the schemes 
with weights (12591 ). (12601) attributed to the problem (12521) . (12531) . where 

A = C + D (271) 

and D, C are specified according to (I267l) - (l269l) or (l267l) . (l268l) . (l270l) . 

Theorem 25. The difference scheme ( 12591) . ( I2g0l) wzY/z (I267l)-(l2g9l). ( 12771) (or 
A267\) , l\268\) , l\270\) , A271\) ) is monotone, and the difference solution satisfies the a 
priori estimate A262\l in L^o (or in Li) under the restriction 

h\v{x±Q.hh,t)\ 

for any t > if a = 1, and ifO<a<l, then this is true under the constraint on 
a time step 

r < \ , (273) 

(1 - (t)7 

with 

7 = max ( ^(k(x + 0.5h) + k(x - 0.5h)) - \(v(x + 0.5h,t) -v(x-0.5h,t))) 
for Ii270\) , and with 

7 = max (^{k{x + 0.5h) + k{x-0.5h)) + ^{v{x + 0.5h,t)-v{x-0.5h,t))) 
x&u) 2h J 

in the case A271\l . 

Proof. Consider the case of the convection-diffusion equation (I243I) - (I245I) (ap- 
proximations (I267l) - (|269l) . (|271l) ) in detail. The problem (I244|) - (I246I) (approxi- 
mations (12671) . (12681) . (12701) ) are investigated is a similar way. 
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To apply Theorems [23] and [241 write explicitly the elements of A. For (1267! )- 
([2691) . ([27T]) ). we have 

- 1 1 1 

ai,i-i - - ^^»~i/2) 

_ 1 , 1 

where fci±i/2 = ± 0.5/i), x G u. 

The condition of nonpositivity of off-diagonal elements (|264l) holds for 

11 11 

^^i"i/2 + ^^^-1/2 > 0' /^^^+i/2 - ^^i+1/2 > 0- (274) 

In this case, diagonal dominance is assured. A spatial computational grid with 
the step from the conditions (12721) satisfies the inequalities (12741) . Restrictions on 
a time step (12611) are reduced to the particular condition (12731 ). Thus, the condi- 
tions of Theorems [23] and [24] hold. This provides the stability and monotonicity 
of the difference solution of the convection-diffusion problem Under the above 
restrictions on the time step. □ 

To overcome restrictions on the spatial grid (12721) . we apply upwind approxi- 
mations for the convective terms. We introduce notation 

v{x, t) = v^{x, t) + t>^(x, t), 
v+{x,t) = ]^{v{x,t) + \v{x,t)\)>Q, 

v-{x,t) = ]^{v{x,t) - \v{xM 



Instead (12691) . we put 

Cw = j^^ i^ + 0.5/i, t){w{x + h,t) — w{x, t)) 

+ j^^^i^ ~ O-S/i, t){w{x, t) — w{x — h, t)). 
For the convective transport in the divergent form, we get 

Cw = ■^('^"(^ + O.S/i, t)w{x + h,t) — v~{x — 0.5h, t)w{x, t)) 

+ ^{v^{x + 0.5/i, t)w{x, t) -v^{x- 0.5/1, t)w{x - h, t)). 



(275) 



(276) 
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Theorem 26. The difference scheme f l259D . ^26U\) with ^267\) . ^269\) . ^271]) . d275D 
( or t\267\l , A268\) , A271\) , l\276\) ) is monotone, and the difference solution satisfies the 
a priori estimate ^2(521) in L^o (or in Li)for any r > if a = 1, and ifO < a < 1, 
then this is true under the constraints on a time step A273\) with 

7 = max (^(k(x + 0.5h) + k(x-0.5h))-]-(v'(x + 0.5h,t)-v^(x-0.5h,t))) 
xglo \h-^ h J 

for ^2751) . and with 

7 = max f-^(A;(x + 0.5/i) + A;(x-0.5/i)) + ^(t;+(x + 0.5/i,t)-t;"(x-0.5/i,t))') 
x&u) h / 

in the case ^276^ . 

In particular, the fully implicit scheme (a = 1) is unconditionally stable and 
monotone. The principal shortcomings of the above schemes are connected with 
the upwind approximations for convective terms (12751) . (12761) ) — these schemes 
indicate the first-order approximation in space. Schemes on the basis of the cen- 
tral difference approximations (12691) . (12701) ) are more accurate — they have the 
second-order spatial approximation. 

5.4. Exponential schemes 

It is convenient to construct monotone schemes by means of transforming the 
original convection-diffusion equation, i.e., by eliminating the convective terns. 
The equation (12431) may be written as 

^ ^ ^k{x)xM^\=f{x,t), {111) 



where 



dt t) dx \ ' dx J 

v{s,t) 



X{x,t) = exp j - j 



k{s) 



ds . (278) 



The equation (12461) is reduced to 

du d f k{x) d{x{x,t)u) 



dt dx \x{x,t) dx J ^ ^ '^^^ 

Further, we can design discretizations in space, i.e., exponential schemes ll7l[33ll. 
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Similarly to (12671 ). for the grid functions satisfying (12681 ). it is possible to put 
in equation (12771 ): 



Aw = —— — -k(x + 0.5h)x(x + 0.5h, t)(w(x + h,t) — w(x, t)) 

k{x — 0.5h)x{x — 0.5h, t){w{x, t) — w{x — h, t)), 



(280) 



h^x{x,t) 

where 

(x~0.5h 

Taking into account that 

(x-0.5h 

with a precision of 0(/i^) we put 

x{x — 0.5h, t) = x{x) exp{9{x, t)h) 

with notation 

Therefore, instead of (|280l) . we can use the following approximation: 

Aw = — -rirkix + 0.5/;.) exp(6'(a;, t)h)(w(x + h,t) — w(x, t)) 

+ -TT^kix — 0.5h) exp(—9(x, t)h)(w(x, t) — wix — h, t)). 
For equation (12791 ). similarly to (|280|) . we put 

Aw = - 7XT~^7rSriA (^(^ + ^' ^)^(^ + h,t)- x{x, t)w{x, t)) 
h'^Xy^ + O.on, t) 

+ - sl^t) t)w{x, t) - x{x - h, t)w{x - h, t)). 



(281) 
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Simplifying this expression, we obtain 



Aw = — -ri^kix + 0.5/;,) exp(— 6'(x + h, t)h)wix + h, t) 
+ T^k^x + 0.5h) exp{9{x, t)h)w{x, t) 

(282) 

+ -r^kix — 0.5/i) exp(— 6'(a;, t)h)w(x, t) 

— —Trkix — 0.5h) exp(^^(x — h, t)h)w(x — h,t). 

Using the above-introduced approximations for the convection-diffusion op- 
erator, we can construct monotone schemes. The primary statement is formulated 
as follows. 

Theorem 27. If on the set of grid functions ^268^ the operator A is defined ac- 
cording to (|2^7]) ( or (|2S2]) ), then the difference scheme ^2591) . Ii260\) . is monotone, 
and the difference solution satisfies the a priori estimate ^262^ in the (or in 
Li)for any r > if a = 1, and ifO < a < 1 then this is true under the constraints 
on a time step ( \273\i with 

7 = max -r7:(k(x + 0.5h) exp(6(x, t)h) + kix — 0.5/i) exp(— ^^(x, t)h)). 

Proof. In the case of (128 1[ for the matrix elements, we have 

= ^(^i+i/2exp(6'i) + A;i_i/2exp(-6'i)), 

tti^i-i = -^ki_i/2exp{-9i), 

aj^j+i = -■^/cj+i/2exp(6'i). 

Checking diagonal dominance by rows and the non-negativity of the off-diagonal 
elements is evident. 

In the case (12821 ). we obtain 

= ^(^i+i/2exp(6'i) + A;i_i/2exp(-6'i)), 

ai,i^i = -^ki^i/2exp{9i-i), 
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exp(-6'i+i). 



In view of the non-negativity of the off-diagonal elements, the condition of diag- 
onal dominance by columns (12571) takes the form 



Thus, the conditions for stability and monotonicity are the same as for schemes 
with the upwind approximations of convective terms (Theorem [261). However, dis- 
cretization in space is of second order as for schemes with the central-difference 
approximations (Theorem [25l). Some complications in evaluating coefficients of 
the difference operator leads to a slight increasing of the computational costs. 

5.5. Multidimensional problems 

Possibilities of constructing second-order monotone schemes for time-dependent 
equations of convection-diffusion are examined on the model 2D problems (|36l)- 
(|I8]) and in the rectangle Vl. 

The convection-diffusion operators in multidimensional problems are repre- 
sented as the sum of the ID convection-diffusion operators. Because of this, in 
constructing monotone schemes for multidimensional problems, we can apply the 
above approximations designed for the ID operators of convection-diffusion. 

Similarly to (I277D . (I278D . rewrite equation ([361) as 



and it is obviously true. 



□ 




(283) 



where now 




(284) 



A similar transformation for (|39l) yields 




(285) 
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For simplicity, we use a uniform grid in each spatial direction. For grids in 
separate directions x^, a = 1,2, we use notation introduced above: 

u = u U duj = coi X 0J2, u = Ui X UJ2. 

After discretization in space of the boundary value problems (l37l) . (l38l) . (12831 
and (I37l), ([38]), (l284l) . we arrive at the problem (l252l) . (12531) . where 

A = Ai + ^2, (286) 

and Aq, a = 1, 2 are ID grid operators of convection-diffusion. On the set of grid 
functions such that 

w(a;,t) = 0, xedu, (287) 
for equation (12831) . similarly to (12811) . we put 

= — -r2k{xi + O.S/ii, X2) exp(6'(a;, t)hi)w{xi + /ii, X2, t) 
"-1 

+ -7-2^(3^1 + 0.5/ii, X2) ex^{6{x, t)hi)w{x, t) 



+ -j^kiyXi — 0.5/ii, X2) exp(— ^^(as, t)hi)w{x, t) 

— -n^k^xi — 0.5/ii, 0:2) exp(— 6'(a;, t)/;,i)w(a;i — /ii, X2, t), 

^42^ = - j^k{xi,X2 + O.5/12) exp(6'(a;, t)h2)w{xi, X2 + h2,t) 
+ -prkixi, X2 + 0.5/12) exp(6'(a;, t)h2)wix, t) 

K 

+ -7-2^(2^15 X2 — 0.5/2,2) exp(— ^^(a3, t)h2)w{x, t) 
1^2 

■^k{xi,X2 - 0.5/12) exp{-9{x, t)h2)w{xi, X2 - h2, t), 



hi 



where 



2k{xj 



(288) 



(289) 
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In the case of (12851) . we have (see (12801 )) 

Aiw = — T2k{xi + O.5/11, X2) exp(— 6'(xi + hi, X2,t)hi)w{xi + hi,X2, t) 



+ T2^(^i + 0.5/ii, X2) exp(^^(x, t)h)w{x, t) 
"■1 



+ —kHxi — 0.5/ii, X2) exp(— 6'(a;, t)hi)w(x, t) 
hi 



— 7-2^(3^1 — 0. 5/2,1, X2) exp{9{xi — hi, X2, , t)h)w{xi — hi, X2, t). 



(290) 



A2W = - j^k{xi,X2 + 0.5/12) exp(-0(xi,X2 + h2,t)hi)w{Xi,X2 + /i2,i) 
+ — j/cfxi, X2 + 0.5/12) exp(^^(£C, t)h)wix, t) 
+ T2^(^i) ^2 — 0.5/12) exp(— ^^(a;, t)hi)w{x, t) 

"■2 

- ^k{xi,X2 - 0.5/12) exp(6'(xi,X2 - h2, ,t)h)w{xi,X2 - h2,t). 
hi 

(291) 

Similarly to Theorem [27l the following statement is proved. 

Theorem 28. If on the set of grid functions ^2^71) the operator A is defined ac- 
cording to ( I2MI) . ( I2MI) . (12191) (or (EMI), ( 12901) . ( 12971) ). difference scheme 
( I259D . ( 12(501) /5 monotone, and the difference solution satisfies the a priori estimate 
A262\) in the L^o (or in Li)for any r > if a = 1, and ifO<a<l, then this is 
true under the constraints on a time step l\272\) with 

7 = max I -r^k{xi + 0.5/ii, X2) exp(^^(a;, t)hi) 

+ ^^(^1 ~ 0.5/ii, X2) exp(— 6'(a;, t)hi) 

+ T2 ^(^1' ^2 + 0.5/12) exp(6'(£c, t)/i2) 

/l2 

+ -hk{xi,X2 - 0.5/12) exp(-6'(a;,t)/i2)|. 
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5.6. Locally one-dimensional schemes 

Computational implementation of the exponential schemes (12591 ). (12601) . (|286l) - 
(12891) and (12591) . (12601) . (12861) . (12871) . (12901) . (1291]) ) involves the inversion of the 
non- self adjoint elliptic grid operators E + arA, where the matrix has strong di- 
agonal dominance either by rows or by columns. To determine the numerical 
solution at a new time level, we can apply iterative methods. Another possibility 
is to use locally one-dimensional schemes, which are based on the splitting (I286|) 
ll44l l27n. Intending to 3D generalizations, we restrict ourselfs to componentwise 
splitting schemes rldblil. 

Rewrite the difference equation (12591 ) as follows: 

yU+l ^ ^ ^^n^ (292) 

where S is the transition operator. For the scheme with weights (12591) . we have 

S={E + aTA)-\E + {(T-l)TA). (293) 

From the stability condition (12601 ). (12921 ). we get 

< 1. (294) 

Monotonicity is ensured by the fact that the matrices {E + arA)^^ and E + {a — 
\)tA are M-matrices. 

Splitting schemes are constructed using transition operators for the individual 
terms in the additive representation (12861) . Let us define 

S^{T) = {E + aTA^)-\E + {a-l)TA^), a = 1,2. (295) 

Instead of (12931 ). we will employ 

S = S,{r)S2{r). (296) 

The stability condition (12941) is true if 

||5J|<1, a = 1,2. (297) 

For the monotonicity of the scheme (12921) . (12961) . it is sufficient to require that the 
individual matrices Sa, a = 1,2 will be M-matrices. For any value of a, only the 
first-order accuracy with respect to r is possible. 
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Numerical implementation of the scheme (12601) . (I292D . (12951) . (12961 ) can be 
conducted using locally one-dimensional schemes with weights, i.e., 

„,n+a/2 _ „,n+(a-l)/2 

^ ^ + A^{ay-^-'' + (1 - ^298) 

= ^l. a = 1,2, 

where, e.g.. 

Theorem 29. Tfon the set of grid functions ( 12^71) operators A^, a = 1,2 
are defined according to A288\) , A289\l (or A280\l , A281\l ), then the locally one- 
dimensional difference scheme A260\) , (\298\) is monotone, and the difference so- 
lution satisfies the a priori estimate A262\) in (or in Li)for any t > Oifa = 1, 
and z/0 < cr < 1, then this is true under the constraints on a time step A273\) with 



7 = max I -i^k{xi + 0.5/ii, X2) exp(6'(a;, t)hi] 
x&uj L hi 



+ -r2k{xi — 0.5/ii, X2) exp(— 6'(a;, t)hi), 
hi 



X2 + 0.5/2,2) exp(6'(a;, t)h2) 

+ ■^k{xi,X2 - 0.5/12) exp(-^(a;,t)/i2)|. 

Proof. Conditions for stability and monotonicity are verified for each individual 
equation (12981) . In particular, for the first equation, we have 

for 

7 = max I -r^f^i^i + 0.5/11,^2) exp(^^(a;, t)hi) 
xeui L hi 

+ -r2^{^i ~ 0.5/ii, 0:2) exp(— 6'(a;, t)hi)\. 
For the second equation, we get 

< ||^n+l/2||^^||^„|| 
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for 



7 = max I -r^^i^ij ^2 + 0.5/12) exp(6'(a;, t)/i2) 

+ T2^(3^i?3;2 - 0.5/12) exp(-6'(a;,t)/;,2)|. 



hi 

Monotonicity of locally one-dimensional schemes under consideration is estab- 
lished in a similar way. □ 

Another classes of splitting schemes can be applied, too. In this regard, we 
highlight the class of additively averaged schemes. 

Instead of a multiplicative representation of the transition operator (12961 ). we 
can employ the additive representation 

5 = ^(5i(2r) + 52(2r)) (299) 

with preserving the first-order approximation in time for the scheme (12921) . 

For the scheme (12921 ). (12951 ). (12991 ). we present another variant of numerical 
implementation. Define the auxiliary functions y^^^, a = 1,2 from 



n+l 

'^"^ +A^{af:-^' + il-a)y2) = 0. (300) 



2r 

For the approximate solution at a new time level, we put 

2/"+' = ^(yr' + 2/2+') + ^^"- (301) 

Conditions of stability and monotonicity for this additively averaged locally 
one-dimensional scheme are formulated in the following theorem. 

Theorem 30. If on the set of grid functions A287\) the operators Aa, a = 1,2 are 
defined according to ( \288\i . A289\i (or ( \290\i . A291\i ). then the additively averaged 
locally one-dimensional difference scheme Ii260\) . A300\) . A301\) . is monotone, and 
the difference solution satisfies the a priori estimate A262\) in Loo (or in Li) for 
any r > if a = 1, and z/O < cr < 1, this is true under the constraints on a time 
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step ({273\) with 



7 = 2 max < —^k(xi + O.bhi, X2) exp(9(x, t)hi) 

x€uj I hf 

+ T^kixi - 0.5hi, X2) exp(-0(x, t)hi), 

K 

■j^k{xi, X2 + 0.5/12) exp(6'(a;, t)h2) 

+ jok{xi,X2 - 0.5/12) exp(-6'(a:;,t)/i2)|. 

Additively average schemes, on the one hand, demonstrate lower accuracy in 
comparison with schemes of componentwise splitting, but on the other hand, they 
are more promising in terms of parallel Computing — the components y^^^, a = 
1, 2 are determined (see (13001) ) independently of each other. 
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